cuda矩阵相乘-02

本博文主要讲解下基于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的，水平有限，有不当之处，请指教，谢谢！