在CUDA中实现更快的矩阵乘法

3

目前,我在cuda c中制作了一个神经网络程序。由于需要操作矩阵乘法,所以我没有使用CUBLAS进行MM。我使用以下代码进行MM。我想知道是否有人有一些建议可以使其更快,这将非常有帮助,因为我需要在学习过程中数百万次使用MM。谢谢。

这是MakeFile:

# cuda root
_CUDA_ROOT_ = /usr/local/cuda

NVCC = nvcc
# include and lib paths
INCLUDES=-I${_CUDA_ROOT_}/include
LIB_PATH=-L${_CUDA_ROOT_}/lib64

# libraries to link against
LIB= -lcudart -lcublas
CU_SRC= main.cu
EXE=$(CU_SRC:.cu=)
#------------------------------
# Choose your gpu arch
SM = sm_35
all: $(EXE)
$(EXE): $(CU_SRC)
        $(NVCC) -arch $(SM) $(CU_SRC) -o $(EXE) $(LIB_PATH) $(LIB)

clean:
        rm -f *.o *.cu_o $(EXE)

这是MM代码:

__global__
void matrixMulti(float* A_d, float* B_d, float* C_d, int m, int k, int n)
{
    __shared__ float ds_A[TILE_WIDTH][TILE_WIDTH];
    __shared__ float ds_B[TILE_WIDTH][TILE_WIDTH];
    int col = blockIdx.x*blockDim.x + threadIdx.x;
    int row = blockIdx.y*blockDim.y + threadIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    float sum = 0;

    for(int t=0; t<(n-1)/TILE_WIDTH+1; t++)
    {
        if(row<m && t*TILE_WIDTH+tx<n)
            ds_A[ty][tx] = A_d[row*n + t*TILE_WIDTH+tx];
        else
            ds_A[ty][tx] = 0.0;
        if(t*TILE_WIDTH+ty<n && col<k)
            ds_B[ty][tx] = B_d[(t*TILE_WIDTH+ty)*k + col];
        else
            ds_B[ty][tx] = 0.0;
        __syncthreads();
        for(int i=0; i<TILE_WIDTH; i++)
            sum += ds_A[ty][i] * ds_B[i][tx];
        __syncthreads();
    }
    if(row<m && col<k)
        C_d[col+row*k] = sum;
}

以下是代码主要部分的示例:

const int TILE_WIDTH = 32;

int main()
{
    int m, k, n;
    m = 10000, k = 10000, n = 10000;
    float *A, *B, *C;
    A = new float[m*n];
    B = new float[n*k];
    C = new float[m*k];
    float *A_d, *B_d, *C_d;
    for (int i=0; i<m*n; i++)
    {
        A[i] = 2;
    }
    for (int i=0; i<n*k; i++)
    {
        B[i] = 3;
    }
    cudaMalloc(&A_d, sizeof(float)*m*n);
    cudaMalloc(&B_d, sizeof(float)*n*k);
    cudaMalloc(&C_d, sizeof(float)*m*k);
    cudaMemcpy(A_d, A, sizeof(float)*m*n, cudaMemcpyHostToDevice);
    cudaMemcpy(B_d, B, sizeof(float)*k*n, cudaMemcpyHostToDevice);
    dim3 dimGrid((k-1)/TILE_WIDTH+1, (m-1)/TILE_WIDTH+1, 1);
    dim3 dimBlock(TILE_WIDTH, TILE_WIDTH, 1);
    matrixMulti<<<dimGrid,dimBlock>>>(A_d, B_d, C_d, m, k, n);
    cudaMemcpy(C, C_d, sizeof(float)*m*k, cudaMemcpyDeviceToHost);
    return 0;
}

3
为什么cuBLAS和BLAS矩阵乘法原语不够好?能否解释一下?您的代码不是很透明。 - Stefano M
这是因为我需要在矩阵乘法代码中添加更多内容,这只是我的代码基础。如果不太清晰,很抱歉。任务是获取C_d = A_d * B_d,其中A_d的维度为mn,B_d的维度为nk。 - user2873565
所使用的算法是瓦片矩阵乘法。 - user2873565
2
让GPU运行快速的一部分是有效的数据组织。如果你强制使用任意数组定义,那么你可能无法有效地使用cublas。相反,你可能想考虑“什么样的数组定义可以让我使用cublas来保存我的权重”。而你在这里展示的代码只不过是“普通的天真瓷砖矩阵乘法”,所以别人很难轻易地理解你的真实意图。cuDNN库是专门为了解决有效地使用类似cublas操作的神经网络权重数组问题而创建的。 - Robert Crovella
1
关于如何偏移指针的后续问题不应该在这里讨论,而应该在一个独立的单独讨论中。然而,A + 100 应该给出指向 A2 的第一个元素的指针... - Jez
显示剩余7条评论
1个回答

10
首先,确保您真的想要这样做。如果没有描述您想要做的操作,很难对此进行评论,但请注意矩阵乘法是一个n立方运算。如果您的操作不具有相同的复杂性,那么使用cuBLAS可能会更好。

为什么呢?cuBLAS可能比您编写的任何内容都要快,并且在GPU架构更新时将更易于维护。像GEMM这样的最佳实现将根据架构而变化,因此您现在针对硬件编写的任何代码都必须重新优化以适应新的硬件。
现在,回答问题。你应该考虑以下几种技术来优化这段代码:

  1. 每个线程计算多个输出值。这减少了共享内存的压力,因为切片数据可以用于多个计算。
  2. 修正共享内存中的冲突。这应该已经被文档详细介绍了。
  3. 向量化共享内存加载和存储。我注意到您正在为sm_35编译。该架构的共享内存库每个时钟周期的带宽为64位。加载单个浮点数只有32位,因此如果不进行向量化,则无法在浮点数上获得完全带宽。您应该查看float2 / float4类型。
  4. 考虑双缓冲。在操作一个共享内存切片时,将数据加载到另一个共享内存切片中。这可以更有效地隐藏全局内存操作的高延迟,并减少同步开销,通常具有更好的性能。但是需要使用两倍的共享内存,因为您需要同时使用两个切片。
对于在GPU上实现矩阵乘法有许多论文,建议您查看它们。从这些论文中获取的详细信息比在SO上提出广泛问题要多得多。
最后... 您确定不想使用cuBLAS吗?我不会指望获得cuBLAS性能的75%,即使那也是一项挑战。

这是这类问题的唯一正确答案。 - talonmies
感谢您的回复。我不使用cuBLAS的主要原因是,例如我在一个数组中拥有所有层(用于神经网络)的权重,并且我需要分离每个层的权重并对每个层进行乘法运算。 在我的当前代码中,这很容易,我只需添加一个for循环和几行代码即可。 但是对于cuBLAS,我找不到一个好的解决方案(如何将大数组的一部分分离出来并将其传递给cuBLAS而不会造成太大的负担)。 请注意,我必须将所有权重放在一个数组中。 - user2873565
你有解决方案吗?如果没有,请给我一些详细的解释,因为它们非常笼统。例如,对于你的第一个解决方案,在当前代码中我应该做什么? - user2873565
所以我有点困惑如何偏移指针。假设我有一个维度为(10,10)的矩阵A1,A2(5,5)和A3(2,2)。现在我把它们都放在数组A中,就像A [A1 A2 A3]一样,这意味着数组A的前100个项目是A1,接下来的25个项目是A2,最后4个项目是A3。 - user2873565
现在,如果我想使用cublasSgemm并访问A2并与某个数组(如B)相乘,我该怎么做?非常感谢。 - user2873565

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接