cuda从入门到精通（六）共享内存和循环分块实现CUDA矩阵乘

本文系转载，出处：https://mp.weixin.qq.com/s/1w1WFPoUEvVECsurqmvJDw
在CUDA编程中，共享内存和循环分块（loop tiling）是两种常见的优化策略，它们可以帮助我们提高矩阵乘法的性能。
共享内存（Shared Memory）：在GPU中，每个线程块（block）都有自己的共享内存。与全局内存相比，共享内存的访问速度更快，但容量较小。因此，如果可能的话，我们应该尽量将数据存储在共享内存中，以减少全局内存访问的延迟。
对于矩阵乘法，我们可以使用共享内存来存储子矩阵的部分结果。每个线程块可以负责计算一个子矩阵的结果，并将结果存储在共享内存中。然后，我们可以使用另一个线程块来将这些子矩阵的结果相加，得到最终的矩阵乘法结果。

循环分块（LoopTiling）：循环分块是将大的循环分解为一系列小的循环，以减少内存访问的冲突和提高内存访问的局部性。在矩阵乘法中，我们可以将大的矩阵分解为一系列小的子矩阵，并分别对每个子矩阵进行乘法运算。

例如，假设我们有一个N×N的矩阵乘法，我们可以将其分解为多个(N/t)×(N/t)的子矩阵乘法，其中t是分块的大小。然后，我们可以使用多个线程块并行计算这些子矩阵的结果，最后将结果相加得到最终的矩阵乘法结果。

下面是一个简单的CUDA代码示例，演示了如何使用共享内存和循环分块来优化矩阵乘法：

__global__ void matMulShared(float* A, float* B, float* C, int N) {
    // 线程块的索引
    int bx = blockIdx.x;
    int by = blockIdx.y;
    // 线程在线程块中的索引
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    // 计算子矩阵的起始位置
    int startRow = N * by;
    int startCol = N * bx;
    // 定义共享内存
    __shared__ float As[tileSize][tileSize];
    __shared__ float Bs[tileSize][tileSize];
    float Csub = 0;
    // 循环分块
    for (int i = startRow; i < startRow + tileSize && i < N; i += tileSize) {
        for (int j = startCol; j < startCol + tileSize && j < N; j += tileSize) {
            // 将子矩阵A和B的数据加载到共享内存中
            for (int m = 0; m < tileSize; m++) {
                As[ty][m] = A[i + m][tx + ty];
                Bs[m][tx] = B[startCol + m][j + tx];
            }
            // 同步线程块中的线程，确保所有线程都加载完数据后再进行计算
            __syncthreads();
            // 计算子矩阵的结果
            for (int m = 0; m < tileSize; m++) {
                Csub += As[ty][m] * Bs[m][tx];
            }
            // 同步线程块中的线程，确保所有线程都计算完结果后再进行下一轮循环
            __syncthreads();
        }
    }
    // 将子矩阵的结果写回全局内存
    int c = startRow * N + startCol + tx + ty;
    if (c < N * N) {
        C[c] = Csub;
    }
}

在上面的代码中，我们使用了tileSize作为分块的大小。As和Bs是两个共享内存数组，用于存储子矩阵A和B的数据。在每个循环迭代中，我们首先将子矩阵A和B的数据加载到共享内存中，然后计算子矩阵的结果，并将结果写回全局内存。我们使用__syncthreads()函数来同步线程块中的线程，确保所有线程都完成了相应的操作后再进行下一轮循环。
请注意，上面的代码只是一个简单的示例，实际上还有很多其他的优化策略和技术可以用来提高矩阵乘法的性能。例如，我们可以使用更复杂的内存访问模式来减少内存访问的冲突，或者使用更高效的算法来计算子矩阵的结果。