CUDA中的非方阵矩阵乘法

3

我需要完成一个非方阵矩阵乘法程序的某些部分,这与我的GPU编程课有关。具体来说,是内核函数和初始化线程块和内核网格的维度。

我参考了CUDA C编程指南中的矩阵乘法代码,但是我修改了我的代码,仅使用给定的参数(因为我们不允许更改参数),而不是像他们一样使用结构体。我们提供了3个矩阵A、B和C,以及它们的维度-m x k、k x n和m x n。在使用结构体时,A.height对应维度m,B.width对应维度n等。

我遇到了几个问题,首先是我的程序未通过包含的测试,该测试验证了乘积矩阵C的正确性。我认为我的矩阵乘法代码存在问题,问题可能出现在我对结构体代码进行修改。

#include <stdio.h>
__global__ void mysgemm(int m, int n, int k, const float *A, const float *B,
        float* C) {

    /********************************************************************
     *
     * Compute C = A x B
     *   where A is a (m x k) matrix
     *   where B is a (k x n) matrix
     *   where C is a (m x n) matrix
     *
     ********************************************************************/

    // INSERT KERNEL CODE HERE
    // Each thread computes one element of C
    // by accumulating results into Cvalue
    float Cvalue = 0;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    for (int e = 0; e < k; ++e){
        Cvalue += (A[row * k + e]) * (B[e * n + col]);
    }
    C[row * n + col] = Cvalue;
}

我另一个问题,我甚至不太确定的是,涉及到初始化线程块和内核网格维度的代码。

// Initialize thread block and kernel grid dimensions ---------------------
    const unsigned int BLOCK_SIZE = 16; // Use 16x16 thread blocks
//INSERT CODE HERE
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid(n / dimBlock.x, m / dimBlock.y);
// Invoke CUDA kernel -----------------------------------------------------
//INSERT CODE HERE
    mysgemm<<<dimGrid, dimBlock>>>(m, n, k, A, B, C);

我理解dimBlock,但我不理解dimGrid,也不知道应该使用什么参数。当我按照原样运行代码时,如果我传递给它的矩阵没有2的幂次维度,内核甚至都不会启动。而如果我使用2的幂,则测试仍然失败。
我很抱歉如果我说得过多了。这是我的第一篇文章,我希望尽可能提供更多的细节。希望有人可以帮忙解决这些问题。

1
有许多关于cuda矩阵乘法的问题,几乎所有可能的变体都被考虑到了。例如,像这个问题一样。也许你应该查看一些已经提出的问题,以获取想法/提示/线索。 - Robert Crovella
2个回答

5
下面的内核是我在CUDA: Tiled matrix-matrix multiplication with shared memory and matrix size which is non-multiple of the block size中发布的内核的变体,不同之处在于它不使用共享内存。
__global__ void MatMulNoShared(float* A, float* B, float* C, int ARows, int ACols, int BRows, int BCols, int CRows, int CCols) {

    float CValue = 0;

    int Row = blockIdx.y*TILE_DIM + threadIdx.y;
    int Col = blockIdx.x*TILE_DIM + threadIdx.x;

    for (int k = 0; k < (TILE_DIM + ACols - 1)/TILE_DIM; k++) {

        for (int n = 0; n < TILE_DIM; ++n) 
            if ((k*TILE_DIM + n < ACols && Row < ARows) && (k*TILE_DIM + n < BRows && Col < BCols))
                CValue += A[Row*ACols + k*TILE_DIM + n] * B[(k*TILE_DIM + n)*BCols + Col];

    }

    if (Row < CRows && Col < CCols) C[((blockIdx.y * blockDim.y + threadIdx.y)*CCols)+(blockIdx.x*blockDim.x)+threadIdx.x]=CValue;
}

内核中的两个if语句是Eric在答案中提到的if语句。
为了方便起见,我在下面发布完整代码:
#include <stdio.h>
#include <math.h>
#include <conio.h>

#define TILE_DIM 16                     // Tile dimension
#define DIMX 373                            
#define DIMY 242
#define DIMZ 533

__global__ void MatMulNoShared(float* A, float* B, float* C, int ARows, int ACols, int BRows, int BCols, int CRows, int CCols) {

    float CValue = 0;

    int Row = blockIdx.y*TILE_DIM + threadIdx.y;
    int Col = blockIdx.x*TILE_DIM + threadIdx.x;

    for (int k = 0; k < (TILE_DIM + ACols - 1)/TILE_DIM; k++) {

        for (int n = 0; n < TILE_DIM; ++n) 
            if ((k*TILE_DIM + n < ACols && Row < ARows) && (k*TILE_DIM + n < BRows && Col < BCols))
                CValue += A[Row*ACols + k*TILE_DIM + n] * B[(k*TILE_DIM + n)*BCols + Col];

    }

    if (Row < CRows && Col < CCols) C[((blockIdx.y * blockDim.y + threadIdx.y)*CCols)+(blockIdx.x*blockDim.x)+threadIdx.x]=CValue;
}

int main() {

    int CCols = DIMZ, CRows=DIMX, ACols=DIMY, ARows=DIMX, BCols=DIMZ, BRows=DIMY;

    dim3 dimBlock(TILE_DIM, TILE_DIM, 1);
    dim3 dimGrid;

    dimGrid.x = (CCols + dimBlock.x - 1)/dimBlock.x;
    dimGrid.y = (CRows + dimBlock.y - 1)/dimBlock.y;

    float *deviceA, *deviceB, *deviceC;

    float* hostA    = (float*)malloc(DIMX*DIMY*sizeof(float));
    float* hostB    = (float*)malloc(DIMY*DIMZ*sizeof(float));
    float* hostC    = (float*)malloc(DIMX*DIMZ*sizeof(float));
    float* hostCp   = (float*)malloc(DIMX*DIMZ*sizeof(float));

    for (int x = 0; x<DIMX; x++)
        for (int y = 0; y<DIMY; y++) {
            hostA[x*DIMY+y] = rand()/(float)RAND_MAX;
            hostB[x*DIMY+y] = rand()/(float)RAND_MAX;
        }

    cudaMalloc((void **)&deviceA, DIMX*DIMY*sizeof(float));
    cudaMalloc((void **)&deviceB, DIMY*DIMZ*sizeof(float));
    cudaMalloc((void **)&deviceC, DIMX*DIMZ*sizeof(float));

    cudaMemcpy(deviceA, hostA, DIMX*DIMY*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(deviceB, hostB, DIMY*DIMZ*sizeof(float), cudaMemcpyHostToDevice);

    MatMulNoShared<<<dimGrid , dimBlock>>>(deviceA , deviceB , deviceC , ARows , ACols, BRows ,BCols , CRows , CCols);

    cudaMemcpy(hostC, deviceC, DIMX*DIMZ*sizeof(float), cudaMemcpyDeviceToHost);

    return 0;
}

请注意这两个指令。
    dimGrid.x = (CCols + dimBlock.x - 1)/dimBlock.x;
    dimGrid.y = (CRows + dimBlock.y - 1)/dimBlock.y;

确保矩阵完全平铺覆盖,如Eric的答案中所提到的第1点。


1

你目前的代码只适用于m和n是16的倍数,也就是你的块大小。

现在有两件事可以做,让它适用于任意大小。

  1. 将网格大小调整为足以覆盖整个矩阵C。不要像您所做的那样使用n/blockdim.x的floor值,可以使用该值的ceil值:

     (n+blockdim.x-1)/blockdim.x
    
  2. 完成步骤1后,由于进行了ceiling操作,您正在乘法的矩阵会稍微变大。您可以通过在内核中添加if子句来限制乘法到结果矩阵C的确切大小。

请参考CUDA文档获取更多详细信息,特别是编程指南。

http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html


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