CUDA加速计算矩阵乘法&进阶玩法(共享内存)

---

一、基础版矩阵乘法

计算下图矩阵P中的一个元素，需要读取矩阵M中对应行的每个元素，矩阵N中对应列的每个元素。单个线程，对读取获得的MN中对应行列数据，进行乘积和操作，从而获得矩阵P中这个元素的结果（简单的矩阵乘法原理）。而GPU并行计算相比较CPU，它可以分配M*N个线程，每个线程负责计算P中的一个元素。

下面为一段利用CUDA计算矩阵乘法的简单示例（注意下面的参数MNK和上图的含义不同，不要误解）

计算矩阵乘法：C = A * B，矩阵A的维度为M*K，矩阵B的维度为K*N

#define M 512
#define K 512
#define N 512

void initial(float *array, int size)
{
	for (int i = 0; i < size; i++)
	{
		array[i] = (float)(rand() % 10 + 1);
	}
}

//核函数（传入显存ABC以及维度信息MNK）
__global__ void multiplicateMatrix(float *array_A, float *array_B, float *array_C, int M_p, int K_p, int N_p)
{
	//这里我们划分的lblock和grid是二维的，分别计算线程的二维索引（x方向和y方向的索引）
	int ix = threadIdx.x + blockDim.x*blockIdx.x;//row number，
	int iy = threadIdx.y + blockDim.y*blockIdx.y;//col number

	if (ix < N_p && iy < M_p)  //筛选线程，每个线程计算C中的一个元素，线程的xy索引与C的元素位置索引对应
	{
		float sum = 0;
		for (int k = 0; k < K_p; k++) //C中的某个元素为A中对应行和B中对应列向量的乘积和。
		{
			sum += array_A[iy*K_p + k] * array_B[k*N_p + ix];
		}
		array_C[iy*N_p + ix] = sum;
	}
}

//主函数
int main(int argc, char **argv)
{
	int Axy = M * K;
	int Bxy = K * N;
	int Cxy = M * N;

	float *h_A, *h_B, *hostRef, *deviceRef;
	
	//在CPU上分配内存
	h_A = (float*)malloc(Axy * sizeof(float)); 
	h_B = (float*)malloc(Bxy * sizeof(float));
	h_C = (float*)malloc(Cxy * sizeof(float));

	initial(h_A, Axy);
	initial(h_B, Bxy);
	
	//在GPU上分配显存
	float *d_A, *d_B, *d_C;
	cudaMalloc((void**)&d_A, Axy * sizeof(float));
	cudaMalloc((void**)&d_B, Bxy * sizeof(float));
	cudaMalloc((void**)&d_C, Cxy * sizeof(float));
	
	//将CPU上初始化的a b值拷贝到GPU上
	cudaMemcpy(d_A, h_A, Axy * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_B, h_B, Bxy * sizeof(float), cudaMemcpyHostToDevice);

	//划分GPU的block和Grid
    int dimx = 2;
    int dimy = 2;
    dim3 block(dimx, dimy);
    dim3 grid((M + block.x - 1) / block.x, (N + block.y - 1) / block.y);
    
    //调用核函数
    multiplicateMatrix<<<grid,block>>> (d_A, d_B, d_C, M, K, N);
    
	//将GPU上计算结果拷贝回CPU
    cudaMemcpy(h_C, d_C, Cxy * sizeof(float), cudaMemcpyDeviceToHost);

	//释放GPU显存资源
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    
	//释放CPU内存资源
    free(h_A);
    free(h_B);
    free(h_C);

return (0);
}

二、为什么可以利用共享内存加速矩阵乘法

一般我们将数据发送到GPU后，默认保存到全局内存，而全局内存的读写速度特别慢，这个时候我们将数据从全局内存放到线程块的共享内存中，计算过程中，读取访问速度更快的共享内存，将会大大减少数据访问耗时，提高程序速度。

1.CUDA内存读写速度比较

下列几种内存的架构参见下图：

• 线程寄存器（~1周期）
• Block共享内存（~5周期）
• Grid全局内存（~500周期）
• Grid常量内存（~5周期）
  

2.申请共享内存

前面对比了共享内存和全局内存的访问速度，为了进一步提高访存速度，可以把全局内存一部分数据拷贝到共享内存中（由于共享内存的大小有限，大概只有几十K，所以只能分多次拷贝数据）

申请共享内存的方式分为静态申请和动态申请
申请共享内存关键字：__ shared __
块内共享内存同步：__syncthreads()函数（块内不同线程之间同步）

• 静态申请

__global__ void staticFun(int* d, int n)
{
	__shared__ int s[64];  //静态申请，需要指定申请内存的大小
	int t = treadIdx.x；
	s[t] = d[t];  //将全局内存数据拷贝到申请的共享内存中，之后利用共享内存中的数据参与运算将会比调
	//用全局内存数据参与运算快(由于共享内存有限，不能全部拷贝到共享内存，这其中就涉及到分批拷贝问题了）
	__syncthreads();//需要等所有线程块都拷贝完成后再进行计算
}

staticFun<<1,n>>(d, n);

• 动态申请

__global__ void dynamicFun(int *d, int n)
{
	extern __shared__ int s[]; //动态申请，不需要指定大小，需要加上extern关键字
	int t = threadIdx.x;
	s[t] = d[t];
	__syncthreads();
}

dynamicFun<<1, n, n*sizeof(int)>>(d, n); //动态申请需要在外部指定共享内存大小

三、改进版矩阵乘法（利用共享内存）

下面的文字很重要，务必结合图进行理解

由于每个Block中，共享内存大小有限，我们可以把矩阵分块计算。
如下图，矩阵Pd中的黄色块，可以由矩阵Md和矩阵Nd中的红框矩阵相乘获得，而红框矩阵的相乘结果，等效于Md中蓝色块与Nd中蓝色块矩阵相乘+Md中橙色块与Nd中橙色块矩阵相乘+…（红框中剩余的分块）。
那么我们计算黄色块，就可以拆分成计算蓝色块矩阵相乘，加上橙色块矩阵相乘，加上其他分块。每次计算，都先将Md和Nd中的分块从全局内存拷贝到共享内存中去。
下面的代码就是依据这个思想进行优化的。

下面直接放源码，读懂上面的图文和下面的代码，可以对共享内存的用法有初步理解。
note:CUDA block的划分即为上图每个分块小矩形的边长（有时为了简化代码我们可以直接划分成正方形）
int dimx = 32;
int dimy = 16;
dim3 block(dimx, dimy);

计算矩阵乘法： C = A * B

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

#define M 512
#define K 512
#define N 512

void initial(float *array, int size)
{
	for (int i = 0; i < size; i++)
	{
		array[i] = (float)(rand() % 10 + 1);
	}
}

//核函数（静态共享内存版）
__global__ void matrixMultiplyShared(float *A, float *B, float *C,
int numARows, int numAColumns, int numBRows, int numBColumns, int numCRows, int numCColumns)
{
	//分配共享内存
	__shared__ float sharedM[blockDim.y][blockDim.x];
	__shared__ float sharedN[blockDim.x][blockDim.y];

	int bx = blockIdx.x;
	int by = blockIdx.y;
	int tx = threadIdx.x;
	int ty = threadIdx.y;

	int row = by * blockDim.y + ty;
	int col = bx * blockDim.x + tx;

	float Csub = 0.0;
	
	//核心：下面将保存在全局内存中的矩阵M&N分块存放到共享内存中
	for (int i = 0; i < (int)(ceil((float)numAColumns / blockDim.x)); i++)//如上图，将一个红框矩形分成多个正方形
	{
		if (i*blockDim.x + tx < numAColumns && row < numARows)//分割M矩阵，边界确定方式结合上图蓝色正方形内数据的位置理解
			sharedM[ty][tx] = A[row*numAColumns + i * blockDim.x + tx];
		else
			sharedM[ty][tx] = 0.0;

		if (i*blockDim.y + ty < numBRows && col < numBColumns)//分割N矩阵
			sharedN[ty][tx] = B[(i*blockDim.y + ty)*numBColumns + col];
		else
			sharedN[ty][tx] = 0.0;
		__syncthreads();//同一线程块中所有线程必须到达运行 __synctrheads()之后才可以做其余操作
		//此操作可以防止当只有部分数据拷贝到共享内存后就提前进行下列计算。

		for (int j = 0; j < blockDim.x; j++)//分块后的矩阵相乘
			Csub += sharedM[ty][j] * sharedN[j][tx];
		__syncthreads();
	}

	if (row < numCRows && col < numCColumns)//将计算后的矩阵块放到结果矩阵C中
		C[row*numCColumns + col] = Csub;
}

//主函数（基本都是常规操作了，和普通版乘法差别不大）
int main(int argc, char **argv)
{
	int Axy = M * K;
	int Bxy = K * N;
	int Cxy = M * N;

	float *h_A, *h_B, *deviceRef;
	h_A = (float*)malloc(Axy * sizeof(float));
	h_B = (float*)malloc(Bxy * sizeof(float));

	deviceRef = (float*)malloc(Cxy * sizeof(float));

	initial(h_A, Axy);
	initial(h_B, Bxy);
	
	float *d_A, *d_B, *d_C;
	cudaMalloc((void**)&d_A, Axy * sizeof(float));
	cudaMalloc((void**)&d_B, Bxy * sizeof(float));
	cudaMalloc((void**)&d_C, Cxy * sizeof(float));

	cudaMemcpy(d_A, h_A, Axy * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_B, h_B, Bxy * sizeof(float), cudaMemcpyHostToDevice);
	
    int dimx = 32;
    int dimy = 16;
    dim3 block(dimx, dimy);
    dim3 grid((M + block.x - 1) / block.x, (N + block.y - 1) / block.y);
	
	matrixMultiplyShared << < grid, block >> > (d_A, d_B, d_C, M, K, K, N, M, N);
	cudaMemcpy(deviceRef, d_C, Cxy * sizeof(float), cudaMemcpyDeviceToHost);
	
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    
    free(h_A);
    free(h_B);
    free(h_C);
}