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(