您的内核中至少有 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;
这应该允许您在以下条件下获得基本正确性:
- 方阵
- 矩阵尺寸可被平铺尺寸整除
修正方阵限制并不难。修复 tile 维度的维度限制需要对内核进行相当大的更改,以便:
- 不处理超出范围的元素
- 使用适合“边界”区域的值正确填充您的共享内存区域
由于您的代码不理解任何这些,我不确定您是否在询问它并选择不专门解决这些问题。
我能够将您的代码的以下改编作为一个基本示例:(请注意,为了减少代码大小,我已经放弃了通常的CUDA 错误检查。请不要将此用作代表性示例良好的编码。进行适当的错误检查。我的回答不是解释良好的 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;
}