矩阵乘法

声明矩阵结构体如下：

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补足。





附：分块矩阵乘法