共享内存矩阵乘法内核

Question

我试图实现基于CUDA C编程指南中概述的基于共享内存的矩阵乘法内核。以下为内核：虽然下面是调用内核

__global__ void matrixMultiplyShared(float * A, float * B, float * C, 
        int ARows, int AColumns, 
        int BRows, int BColumns, 
        int CRows, int CColumns) { 
    float * CSub = &C[CColumns * 16 * blockIdx.y + 16 * blockIdx.x]; 
    float CValue = 0; 
for (int k = 0; k < (AColumns/16); ++k) { 
     float * ASub = &A[AColumns * 16 * blockIdx.y + 16 * k]; 
     float * BSub = &B[AColumns*16*k + 16*blockIdx.y]; 
     __shared__ float As[16][16]; 
     __shared__ float Bs[16][16]; 
     As[threadIdx.y][threadIdx.x] = ASub[threadIdx.y*AColumns+threadIdx.x]; 
     Bs[threadIdx.y][threadIdx.x] = BSub[threadIdx.y*AColumns+threadIdx.x]; 
     __syncthreads(); 
     for (int n = 0; n < 16; ++n) 
     CValue += As[threadIdx.y][n] * Bs[n][threadIdx.x]; 
     __syncthreads(); 
    } 
    CSub[threadIdx.x*CColumns+threadIdx.y]=CValue; 
} 

：

dim3 dimBlock(16, 16, 1); 
dim3 dimGrid; 
dimGrid.x = (CColumns + dimBlock.x - 1)/dimBlock.x; 
dimGrid.y = (CRows + dimBlock.y - 1)/dimBlock.y; 
matrixMultiplyShared<<<dimGrid , dimBlock>>>(deviceA , deviceB , deviceC , ARows , AColumns, BRows ,BColumns , CRows , CColumns); 

不幸的是，这似乎产生不正确的结果。

任何援助/解释将不胜感激。

Answer 1

你的内核至少有2个基本错误，两者都相对简单。当你有这样的：

 float * BSub = &B[AColumns*16*k + 16*blockIdx.y]; 

你应该使用这样的：

 float * BSub = &B[AColumns*16*k + 16*blockIdx.x]; 

而且，你有这样的：

CSub[threadIdx.x*CColumns+threadIdx.y]=CValue; 

你应该使用这样的：

CSub[threadIdx.y*CColumns+threadIdx.x]=CValue; 

这应该让你获得基本的正确性在以下条件下：

1. 方阵
2. 矩阵尺寸由瓦片尺寸

固定方阵限制整除并不困难。固定尺寸的限制对瓷砖尺寸涉及到内核相当大的变化，依次是：

1. 不处理超出范围的元素
2. 正确地是在适当的值填充您的共享内存区域中的“边界“地区

由于您的代码不能理解任何这些，所以我不确定您是否在问这个问题，并选择不专门解决这些问题。

我能得到你的代码的工作作为一个基本的例子如下调整： （请注意，减少代码大小看的利益，我已经免除了通常CUDA error checking请不要使用这个作为。良好的编码的一个代表性的例子。做适当的错误检查。我的回答的一点是不要解释好CUDA错误检查，但表现出的算法正确的例子。）

#include <stdio.h> 
#include <math.h> 
#define TILE_DIM 16 
#define DIMX 256 
#define DIMY 256 
#define RES 0.1 

__global__ void matrixMultiplyShared(float * A, float * B, float * C, 
        int ARows, int AColumns, 
        int BRows, int BColumns, 
        int CRows, int CColumns) { 
    float CValue = 0; 
    if (((blockIdx.y * blockDim.y + threadIdx.y)< CRows) && ((blockIdx.x * blockDim.x + threadIdx.x) < CColumns)) { 
     for (int k = 0; k < (AColumns/TILE_DIM); ++k) { 
     float * ASub = &A[AColumns * TILE_DIM * blockIdx.y + TILE_DIM * k]; 
     float * BSub = &B[AColumns*TILE_DIM*k + TILE_DIM*blockIdx.x]; 
     __shared__ float As[TILE_DIM][TILE_DIM]; 
     __shared__ float Bs[TILE_DIM][TILE_DIM]; 
     As[threadIdx.y][threadIdx.x] = ASub[threadIdx.y*AColumns+threadIdx.x]; 
     Bs[threadIdx.y][threadIdx.x] = BSub[threadIdx.y*AColumns+threadIdx.x]; 
     __syncthreads(); 
     for (int n = 0; n < TILE_DIM; ++n) 
     CValue += As[threadIdx.y][n] * Bs[n][threadIdx.x]; 
     __syncthreads(); 
     } 
     C[((blockIdx.y * blockDim.y + threadIdx.y)*CColumns)+(blockIdx.x*blockDim.x)+threadIdx.x]=CValue; 
    } 
} 


void matrixMultiplyCPU(float * A, float * B, float * C, 
        int ARows, int AColumns, 
        int BRows, int BColumns, 
        int CRows, int CColumns) { 
    for (int i = 0; i<ARows; i++) 
    for (int j=0; j<BColumns; j++){ 
     float Ctemp = 0.0; 
     for (int k=0; k<AColumns; k++) 
     Ctemp += A[i*AColumns + k] * B[k*BColumns+j]; 
     C[i*CColumns+j] = Ctemp; 
     } 

} 
int main(){ 
int CColumns = DIMY, CRows=DIMX, AColumns=DIMY, ARows=DIMX, BColumns=DIMY, BRows=DIMX; 
dim3 dimBlock(TILE_DIM, TILE_DIM, 1); 
dim3 dimGrid; 
dimGrid.x = (CColumns + dimBlock.x - 1)/dimBlock.x; 
dimGrid.y = (CRows + dimBlock.y - 1)/dimBlock.y; 
float *deviceA, *deviceB, *deviceC; 
float hostA[DIMY][DIMX]; 
float hostB[DIMY][DIMX]; 
float hostC[DIMY][DIMX]; 
float hostCp[DIMY][DIMX]; 
for (int x = 0; x<DIMX; x++) 
    for (int y = 0; y<DIMY; y++) { 
    hostA[y][x] = rand()/(float)RAND_MAX; 
    hostB[y][x] = rand()/(float)RAND_MAX; 
    } 
    cudaMalloc((void **)&deviceA, DIMX*DIMY*sizeof(float)); 
    cudaMalloc((void **)&deviceB, DIMX*DIMY*sizeof(float)); 
    cudaMalloc((void **)&deviceC, DIMX*DIMY*sizeof(float)); 
    cudaMemcpy(deviceA, hostA, DIMX*DIMY*sizeof(float), cudaMemcpyHostToDevice); 
    cudaMemcpy(deviceB, hostB, DIMX*DIMY*sizeof(float), cudaMemcpyHostToDevice); 
    matrixMultiplyShared<<<dimGrid , dimBlock>>>(deviceA , deviceB , deviceC , ARows , AColumns, BRows ,BColumns , CRows , CColumns); 
    cudaMemcpy(hostC, deviceC, DIMX*DIMY*sizeof(float), cudaMemcpyDeviceToHost); 
    matrixMultiplyCPU(&(hostA[0][0]) , &(hostB[0][0]) , &(hostCp[0][0]) , ARows , AColumns, BRows ,BColumns , CRows , CColumns); 

for (int y = 0; y<DIMY; y++) 
    for (int x = 0; x<DIMX; x++) 
    if (fabs(hostCp[y][x] - hostC[y][x]) > RES) 
     { 
     printf("Error at offset y=%d,x=%d, CPU = %f, GPU = %f\n", y, x, hostCp[y][x], hostC[y][x]); 
     return 1; 
     } 
printf("Finished!\n"); 
return 0; 
}