声明矩阵结构体如下:
typedef struct{
int width;
int height;
float* elements;
}Matrix;
串行版本:
void MatMul(const Matrix A,const Matrix B,Matrix C)
{
for(int row=0;row<C.height;row++)
{
for(int col=0;col<C.width;col++)
{
float Cvalue=0;
for(int e=0;e<A.width;e++)
{
Cvalue+=A.elements[row*A.width+e]*B.elements[e*B.width+col];
}
C.elements[row*C.width+col]=Cvalue;
}
}
cuda基本版本:
#define BLOCK_SIZE 16
__global__ void MatMulKernel(Matrix A,Matrix B,Matrix C){
float Cvalue=0;
int row=blockIdx.y*blockDim.y+threadIdx.y;
int col=blockIdx.x*blockDim.x+threadIdx.x;
if(row<C.height&&col<C.width)
{
for(int e=0;e<A.width;++e)
Cvalue+=A.elements[row*A.width+e]*B.elements[e*B.width+col];
C.elements[row*C.width+col]=Cvalue;
}
}
void MatMul(const Matrix A,const Matrix B,Matrix C)
{
Matrix d_A;
d_A.width=A.width;d_A.height=A.height;
size_t size=A.width*A.height*sizeof(float);
cudaMalloc((void**)&d_A.elements,size);
cudaMemcpy(d_A.elements,A.elements,size,
cudaMemcpyHostToDevice);
Matrix d_B;
d_B.width=B.width;d_B.height=B.height;
size=B.width*B.height*sizeof(float);
cudaMalloc((void**)&d_B.elements,size);
cudaMemcpy(d_B.elements,B.elements,size,
cudaMemcpyHostToDevice);
Matrix d_C;
d_C.width=C.width;d_C.height=C.height;
size=C.width*C.height*sizeof(float);
cudaMalloc((void**)&d_C.elements,size);
dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
dim3 dimGrid((B.width-1)/dimBlock.x+1,(A.height-1)/dimBlock.y+1);
MatMulKernel<<<dimGrid,dimBlock>>>(d_A,d_B,d_C);
cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
cuda优化版本
利用shared memory,减少对内存的访问次数,从而进一步加速。
代码注释:
9到11行,为分块矩阵Mds和Nds赋值
12到14行,求分块矩阵Mds和Nds的乘积
8行for循环,保证取完Md的一行和Nd的一列,以便完成结果矩阵的计算。eg:Pd00=Mds00*Nds00+Mds01*Nds10
在cuda计算中,矩阵分块的大小是一样的,不足的元素用0补足。
附:分块矩阵乘法