CUDA/GPU下矩阵乘法的几种实现的C++源码

#include "cutil_inline.h"
#include "cublas.h"

void MatrixMul(const float *A, const float *B, float *C, int Width) {
int i, j, k;
for(i=0; i<Width; i++)
for(j=0; j<Width; j++){
float s=0;
for(k=0; k<Width; k++)
s+=A[i*Width+k]*B[k*Width+j];
C[i*Width+j]=s;
}
}

#define TILE_WIDTH 16
//简单方法
__global__ void MatrixMulKernel1(float* Md, float* Nd, float* Pd, int Width)
{
int x = threadIdx.x+blockIdx.x*blockDim.x;
int y = threadIdx.y+blockIdx.y*blockDim.y;

float Pvalue = 0;
for (int k = 0; k < Width; ++k)
Pvalue+=Md[y * Width + k]*Nd[k * Width + x];
Pd[y*Width + x] = Pvalue;
}

//块方法
__global__ void MatrixMulKernel2(float* Md, float* Nd, float* Pd, int Width)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float Pvalue = 0;

for (int m = 0; m < gridDim.x; ++m) {
__shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
__shared__ float Nds[TILE_WIDTH][TILE_WIDTH];
Mds[ty][tx] = *(Md + (by*blockDim.y + ty) * Width + m*blockDim.x + tx);
Nds[ty][tx] = *(Nd + (m*blockDim.y + ty) * Width + bx*blockDim.x + tx);
__syncthreads();
for (int k = 0; k < blockDim.x; ++k)
Pvalue += Mds[ty][k] * Nds[k][tx];
__syncthreads();
}
Pd[(by*blockDim.y+ty)*Width+bx*blockDim.x+tx] = Pvalue;
}
//块方法, 循环展开
__global__ void MatrixMulKernel3(float* Md, float* Nd, float* Pd, int Width)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float Pvalue = 0;

for (int m = 0; m < gridDim.x; ++m) {
__shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
__shared__ float Nds[TILE_WIDTH][TILE_WIDTH];
Mds[ty][tx] = *(Md + (by*blockDim.y + ty) * Width + m*blockDim.x + tx);
Nds[ty][tx] = *(Nd + (m*blockDim.y + ty) * Width + bx*blockDim.x + tx);
__syncthreads();
Pvalue += Mds[ty][0] * Nds[0][tx] + Mds[ty][1] * Nds[1][tx] +
Mds[ty][2] * Nds[2][tx] + Mds[ty][3] * Nds[3][tx] +
Mds[ty][4] * Nds[4][tx] + Mds[ty][5] * Nds[5][tx] +
Mds[ty][6] * Nds[6][tx] + Mds[ty][7] * Nds[7][tx] +
Mds[ty][8] * Nds[8][tx] + Mds[ty][9] * Nds[9][tx] +
Mds[ty][10] * Nds[10][tx] + Mds[ty][11] * Nds[11][tx] +
Mds[ty][12] * Nds[12][tx] + Mds[ty][13] * Nds[13][tx] +
Mds[ty][14] * Nds[14][tx] + Mds[ty][15] * Nds[15][tx];
__syncthreads();
}
Pd[(by*blockDim.y+ty)*Width+bx*blockDim.x+tx] = Pvalue;
}

#define ThreadColumn 2
//块方法, 线程粒度
__global__ void MatrixMulKernel4(float* Md, float* Nd, float* Pd, int Width)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float Pvalue1 = 0, Pvalue2=0;

for (int m = 0; m < gridDim.x; ++m) {
__shared__ float Mds[TILE_WIDTH][ThreadColumn*TILE_WIDTH];
__shared__ float Nds[ThreadColumn*TILE_WIDTH][ThreadColumn*TILE_WIDTH];
for(int k=0;k<ThreadColumn; ++k)
Mds[ty][tx+k*TILE_WIDTH] = *(Md + (by*blockDim.y + ty) * Width +  ThreadColumn*m*blockDim.x + tx + k*TILE_WIDTH );
for(int h=0;h<ThreadColumn; ++h)
for(int k=0; k<ThreadColumn; ++k)
Nds[ty+h*TILE_WIDTH][tx+k*TILE_WIDTH] = *(Nd + (m*blockDim.y*ThreadColumn + ty + h*TILE_WIDTH) * Width  +
bx*blockDim.x*ThreadColumn + tx + k*TILE_WIDTH);
__syncthreads();
for(int h=0;h<ThreadColumn*TILE_WIDTH;++h)
Pvalue1 +=Mds[ty][h] * Nds[h][tx];
for(int k=0;k<ThreadColumn*TILE_WIDTH;++k)
Pvalue2 +=Mds[ty][k] * Nds[k][tx+TILE_WIDTH];
__syncthreads();
}

Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx] = Pvalue1;
Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx + TILE_WIDTH] = Pvalue2;
}
//块方法, 循环展开, 线程粒度
__global__ void MatrixMulKernel5(float* Md, float* Nd, float* Pd, int Width)
{
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
float Pvalue1 = 0, Pvalue2=0;

for (int m = 0; m < gridDim.x; ++m) {
__shared__ float Mds[TILE_WIDTH][ThreadColumn*TILE_WIDTH];
__shared__ float Nds[ThreadColumn*TILE_WIDTH][ThreadColumn*TILE_WIDTH];
for(int k=0;k<ThreadColumn; ++k)
Mds[ty][tx+k*TILE_WIDTH] = *(Md + (by*blockDim.y + ty) * Width +  ThreadColumn*m*blockDim.x + tx + k*TILE_WIDTH );
for(int h=0;h<ThreadColumn; ++h)
for(int k=0; k<ThreadColumn; ++k)
Nds[ty+h*TILE_WIDTH][tx+k*TILE_WIDTH] = *(Nd + (m*blockDim.y*ThreadColumn + ty + h*TILE_WIDTH) * Width  +
bx*blockDim.x*ThreadColumn + tx + k*TILE_WIDTH);
__syncthreads();
Pvalue1 += Mds[ty][0] * Nds[0][tx] + Mds[ty][1] * Nds[1][tx] +
Mds[ty][2] * Nds[2][tx] + Mds[ty][3] * Nds[3][tx] +
Mds[ty][4] * Nds[4][tx] + Mds[ty][5] * Nds[5][tx] +
Mds[ty][6] * Nds[6][tx] + Mds[ty][7] * Nds[7][tx] +
Mds[ty][8] * Nds[8][tx] + Mds[ty][9] * Nds[9][tx] +
Mds[ty][10] * Nds[10][tx] + Mds[ty][11] * Nds[11][tx] +
Mds[ty][12] * Nds[12][tx] + Mds[ty][13] * Nds[13][tx] +
Mds[ty][14] * Nds[14][tx] + Mds[ty][15] * Nds[15][tx] +
Mds[ty][16] * Nds[16][tx] + Mds[ty][17] * Nds[17][tx] +
Mds[ty][18] * Nds[18][tx] + Mds[ty][19] * Nds[19][tx] +
Mds[ty][20] * Nds[20][tx] + Mds[ty][21] * Nds[21][tx] +
Mds[ty][22] * Nds[22][tx] + Mds[ty][23] * Nds[23][tx] +
Mds[ty][24] * Nds[24][tx] + Mds[ty][25] * Nds[25][tx] +
Mds[ty][26] * Nds[26][tx] + Mds[ty][27] * Nds[27][tx] +
Mds[ty][28] * Nds[28][tx] + Mds[ty][29] * Nds[29][tx] +
Mds[ty][30] * Nds[30][tx] + Mds[ty][31] * Nds[31][tx];

Pvalue2 += Mds[ty][0] * Nds[0][tx+TILE_WIDTH] + Mds[ty][1] * Nds[1][tx+TILE_WIDTH] +
Mds[ty][2] * Nds[2][tx+TILE_WIDTH] + Mds[ty][3] * Nds[3][tx+TILE_WIDTH] +
Mds[ty][4] * Nds[4][tx+TILE_WIDTH] + Mds[ty][5] * Nds[5][tx+TILE_WIDTH] +
Mds[ty][6] * Nds[6][tx+TILE_WIDTH] + Mds[ty][7] * Nds[7][tx+TILE_WIDTH] +
Mds[ty][8] * Nds[8][tx+TILE_WIDTH] + Mds[ty][9] * Nds[9][tx+TILE_WIDTH] +
Mds[ty][10] * Nds[10][tx+TILE_WIDTH] + Mds[ty][11] * Nds[11][tx+TILE_WIDTH] +
Mds[ty][12] * Nds[12][tx+TILE_WIDTH] + Mds[ty][13] * Nds[13][tx+TILE_WIDTH] +
Mds[ty][14] * Nds[14][tx+TILE_WIDTH] + Mds[ty][15] * Nds[15][tx+TILE_WIDTH] +
Mds[ty][16] * Nds[16][tx+TILE_WIDTH] + Mds[ty][17] * Nds[17][tx+TILE_WIDTH] +
Mds[ty][18] * Nds[18][tx+TILE_WIDTH] + Mds[ty][19] * Nds[19][tx+TILE_WIDTH] +
Mds[ty][20] * Nds[20][tx+TILE_WIDTH] + Mds[ty][21] * Nds[21][tx+TILE_WIDTH] +
Mds[ty][22] * Nds[22][tx+TILE_WIDTH] + Mds[ty][23] * Nds[23][tx+TILE_WIDTH] +
Mds[ty][24] * Nds[24][tx+TILE_WIDTH] + Mds[ty][25] * Nds[25][tx+TILE_WIDTH] +
Mds[ty][26] * Nds[26][tx+TILE_WIDTH] + Mds[ty][27] * Nds[27][tx+TILE_WIDTH] +
Mds[ty][28] * Nds[28][tx+TILE_WIDTH] + Mds[ty][29] * Nds[29][tx+TILE_WIDTH] +
Mds[ty][30] * Nds[30][tx+TILE_WIDTH] + Mds[ty][31] * Nds[31][tx+TILE_WIDTH];
__syncthreads();
}

Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx] = Pvalue1;
Pd[(by*blockDim.y + ty) * Width +  bx*ThreadColumn*blockDim.x + tx + TILE_WIDTH] = Pvalue2;
}
int main(int argc, char **argv) {
int Width(0);
bool is_test=false;
if(argc==1){
Width=512;
is_test=false;
}
else if(argc==2){
Width=atoi(argv[1]);
is_test=false;
}
else{
Width=atoi(argv[1]);
is_test=(bool)atoi(argv[2]);
}
float *h_A=(float*)malloc(Width*Width*sizeof(float));
float *h_B=(float*)malloc(Width*Width*sizeof(float));
float *h_C=(float*)malloc(Width*Width*sizeof(float));
float *h_C_ref=(float*)malloc(Width*Width*sizeof(float));
float *d_A, *d_B, *d_C;
float t0, error_norm=0, ref_norm=0;
unsigned int timer1=0;
cutCreateTimer(&timer1);
cutStartTimer(timer1);
for(int i=0; i<Width*Width; i++) {
h_A[i]=rand()/(float)RAND_MAX;
h_B[i]=rand()/(float)RAND_MAX;
}
cublasInit();
cublasAlloc(Width*Width, sizeof(float), (void**)&d_A);
cublasAlloc(Width*Width, sizeof(float), (void**)&d_B);
cublasAlloc(Width*Width, sizeof(float), (void**)&d_C);
cublasSetVector(Width*Width, sizeof(float), h_A, 1, d_A, 1);
cublasSetVector(Width*Width, sizeof(float), h_B, 1, d_B, 1);

if(is_test)
MatrixMul(h_A, h_B, h_C_ref,Width);

dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);
dim3 dimGrid(Width/TILE_WIDTH, Width/TILE_WIDTH);

t0=cutGetTimerValue(timer1);
cudaThreadSynchronize();
MatrixMulKernel1<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);
cudaThreadSynchronize();
float gpu_t1=(cutGetTimerValue(timer1)-t0)/1000.0f;
cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);

if(is_test){
error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
float diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}
t0=cutGetTimerValue(timer1);
cudaThreadSynchronize();
MatrixMulKernel2<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);
cudaThreadSynchronize();
float gpu_t2=(cutGetTimerValue(timer1)-t0)/1000.0f;
cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);
if(is_test){
error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
float diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

t0=cutGetTimerValue(timer1);
cudaThreadSynchronize();
MatrixMulKernel3<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);
cudaThreadSynchronize();
float gpu_t3=(cutGetTimerValue(timer1)-t0)/1000.0f;
cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);
if(is_test){
error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
float diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

dimGrid.x/=ThreadColumn;

t0=cutGetTimerValue(timer1);
cudaThreadSynchronize();
MatrixMulKernel4<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);
cudaThreadSynchronize();
float gpu_t4=(cutGetTimerValue(timer1)-t0)/1000.0f;
cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);
if(is_test){
error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
float diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

t0=cutGetTimerValue(timer1);
cudaThreadSynchronize();
MatrixMulKernel5<<<dimGrid,dimBlock>>>(d_A,d_B,d_C,Width);
cudaThreadSynchronize();
float gpu_t5=(cutGetTimerValue(timer1)-t0)/1000.0f;
cublasGetVector(Width*Width, sizeof(float), d_C, 1, h_C, 1);
if(is_test){
error_norm=0;
ref_norm=0;
for(int i=0; i<Width*Width; i++) {
float diff=h_C_ref[i]-h_C[i];
error_norm+=diff*diff;
ref_norm+=h_C_ref[i]*h_C_ref[i];
}
printf("Test %s/n", (sqrtf(error_norm/ref_norm)<1E-6) ? "PASSED" : "FAILED");
}

printf("矩阵阶数为%4d,", Width);
printf("简单方法: %.6fs(%.3fGflops),", gpu_t1, 1e-9*Width*Width*Width*2/gpu_t1);
printf("块方法: %.6fs(%.3fGflops),", gpu_t2, 1e-9*Width*Width*Width*2/gpu_t2);
printf("块+循环展开方法: %.6fs(%.3fGflops),", gpu_t3, 1e-9*Width*Width*Width*2/gpu_t3);
printf("块+线程粒度2: %.6fs(%.3fGflops),", gpu_t4, 1e-9*Width*Width*Width*2/gpu_t4);
printf("块+循环展开方法+线程粒度2: %.6fs(%.3fGflops)/n", gpu_t5, 1e-9*Width*Width*Width*2/gpu_t5);
}

环境：CUDA toolkit3.2+Windows XP+CUDA SDK中的vs2008模板release编译通过，显卡是GeForce GT240。大家可根据自己的情况进行测试，报告一下结果吧.

不同规模的运行结果:
矩阵阶数为 512,简单方法: 0.033887s(7.922Gflops),块方法: 0.008424s(31.865Gflops),块+循环展开方法: 0.003995s(67.191Gflops),块+线程粒度2: 0.008091s(33.179Gflops),块+循环展开方法+线程粒度2: 0.003399s(78.984Gflops)

矩阵阶数为1024,简单方法: 0.264424s(8.121Gflops),块方法: 0.066201s(32.439Gflops),块+循环展开方法: 0.030827s(69.662Gflops),块+线程粒度2: 0.063264s(33.945Gflops),块+循环展开方法+线程粒度2: 0.026282s(81.710Gflops)

矩阵阶数为2048,简单方法: 2.112513s(8.132Gflops),块方法: 0.527002s(32.599Gflops),块+循环展开方法: 0.244058s(70.392Gflops),块+线程粒度2: 0.505495s(33.986Gflops),块+循环展开方法+线程粒度2: 0.208845s(82.261Gflops)

矩阵阶数为4096,简单方法: 17.705070s(7.763Gflops),块方法: 4.205308s(32.682Gflops),块+循环展开方法: 1.966043s(69.906Gflops),块+线程粒度2: 4.059064s(33.860Gflops),块+循环展开方法+线程粒度2: 1.677771s(81.918Gflops)

测试结果表明在块大小16x16时，块方法和循环展开对速度有很显著的影响，而线程粒度的使用对速度只有很小的提高。而块大小是8x8时的情形没有测试，线程粒度对速度应该有较大影响，家可以做一做。

参考文献: David B. Kirk, Wen-mei W. Hwu. 大规模并行处理器编程实战[M]. 北京: 清华大学出版社, 2010.