CUDA系列三：矩阵相乘

本博文主要讲解下基于cuda的矩阵相乘，cuda特别擅长的就是矩阵乘法，而且也比较容易实现。通过矩阵乘法的实现，可以比较容易理解cuda的核心思想。网上也有很多基于cuda实现的矩阵乘法，但是感觉都不完成，要不就是有错，本文给出的代码都是经过验证可行的，希望能够帮助到大家。

矩阵乘法实现方式一：矩阵乘法的逐点实现方式，具体如下图所示

对应实现代码：

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>


__global__ void MatMul(int *M,int *N,int *P,int width)
{
	int x = threadIdx.x;
	int y = threadIdx.y;
	
	float Pervalue = 0;
	
	float elem1 = 0.0,elem2 = 0.0,value = 0.0;
	for(int i = 0;i < width;i++)
	{
		elem1 = M[y * width + i];//取M矩阵的一行
		elem2 = N[i * width + x];//取N矩阵的一列
		
		value += elem1 * elem2;//求和
	}
	
	P[y * width + x] = value;
}

int main()
{
	const int ND = 30;
	int a[ND][ND],b[ND][ND],c[ND][ND];
	int *M,*N,*P;
	
	int width = ND;
	int NUM = 900;
	dim3 blockSize(ND,ND);
	
	cudaEvent_t start,stop;
	float elapsedTime = 0;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	
	//设备端内存分配
	cudaMalloc((void**)&M,ND * ND * sizeof(int));
	cudaMalloc((void**)&N,ND * ND * sizeof(int));
	cudaMalloc((void**)&P,ND * ND * sizeof(int));
	
	//初始化
	for(int i = 0;i < ND;i++)
	{
		for(int j = 0;j < ND;j++)
		{
			a[i][j] = 2;
			b[i][j] = 3;
		}
	}
	
	int Size = ND * ND;
	//数据拷贝，主机到设备
	cudaMemcpy(M,a,Size * sizeof(int),cudaMemcpyHostToDevice);
	cudaMemcpy(N,b,Size * sizeof(int),cudaMemcpyHostToDevice);
	
	cudaEventRecord(start,0);
	MatMul<<<1,blockSize>>>(M,N,P,width);//调用核函数
	cudaThreadSynchronize();
	cudaEventRecord(stop,0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime,start,stop);
	
	cudaMemcpy(c,P,Size * sizeof(int),cudaMemcpyDeviceToHost);
	
	printf("c0 = %d \n",c[0][0]);
	
	//释放设备内存
	cudaFree(M);
	cudaFree(N);
	cudaFree(P);
	
	return 0;
}

运行结果：

矩阵相乘实现方式二：矩阵乘法分块实现，具体如下图所示

具体代码实现：

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>


#define TILE_WIDTH 10

//核函数的具体实现
__global__ void matmul(int *M,int *N,int *P,int width)
{
	__shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
	__shared__ float Nds[TILE_WIDTH][TILE_WIDTH];
	
	int bx = blockIdx.x;
	int by = blockIdx.y;
	int tx = threadIdx.x;
	int ty = threadIdx.y;
	
	int Col = bx * TILE_WIDTH + tx;
	int Row = by * TILE_WIDTH + ty;
	
	int Pervalue = 0;
	
	for(int i = 0;i < width / TILE_WIDTH;i++)  //有多少个TILE_WIDTH，每个循环计算一个块的大小
	{
		Mds[ty][tx] = M[Row * width + (i * TILE_WIDTH + tx)];
		Nds[ty][tx] = N[Col + (i * TILE_WIDTH + ty) * width];
		__syncthreads();
		
		
		for(int k = 0;k < TILE_WIDTH;k++) //TILE_WIDTH相乘
			Pervalue += Mds[ty][k] * Nds[k][tx];
		__syncthreads();
	}
	
	P[Row * width + Col] = Pervalue;
}


int main()
{
	const int Nd = 30;
	int Size = Nd * Nd;
	int *M,*N,*P;
	int width = Nd / 3;
	
	int a[Nd][Nd];
	int b[Nd][Nd];
	int c[Nd][Nd];
	
	//线程块以及线程的划分
	dim3 gridSize(Nd / width,Nd / width);
	dim3 blockSize(width,width);
	
	cudaEvent_t start,stop;
	float elapsedTime;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	
	//设备内存分配
	cudaMalloc((void**)&M,Size * sizeof(int));
	cudaMalloc((void**)&N,Size * sizeof(int));
	cudaMalloc((void**)&P,Size * sizeof(int));
	
	//初始化
	for(int i = 0;i < Nd;i++)
	{
		for(int j = 0;j < Nd;j++)
		{
			a[i][j] = 2;
			b[i][j] = 3;
		}
	}
	
	//数据拷贝，主机到设备
	cudaMemcpy(M,a,Size * sizeof(int),cudaMemcpyHostToDevice);
	cudaMemcpy(N,b,Size * sizeof(int),cudaMemcpyHostToDevice);
	
	cudaEventRecord(start,0);
	matmul<<<gridSize,blockSize>>>(M,N,P,Nd); //调用核函数
	cudaThreadSynchronize();
	cudaEventRecord(stop,0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime,start,stop);
	
	
	cudaMemcpy(c,P,Size * sizeof(int),cudaMemcpyDeviceToHost);
	printf("c0 = %d\n",c[0][0]);
	
	
	cudaFree(M);
	cudaFree(N);
	cudaFree(P);
	
	return 0;
}

运行结果：

本文也参考了网上的一些资料，主要是做了一定的修改以及程序的完备，图片就直接网上copy的，水平有限，有不当之处，请指教，谢谢！