【CUDA】三、存储器层级结构及共享存储器

一、存储器简介

        对于CUDA线程，主要设置的存储器有：局部存储器、共享存储器、全局存储器、只读存储器。

        局部存储器：作用域在每个线程的内部，每个线程都有一个自己私有的存储器；

        共享存储器：作用域在线程块内，相同块内的线程可访问该存储中的数据，已达到数据共享；

        全局存储器：比共享存储器的作用域更大，所有块中的线程都可用访问该存储器中的数据

        只读存储器：包括常量存储器和纹理存储器

其中，全局存储器、只读存储器（常量存储器、文理存储器）的生命周期同当前应用程序的当前内核发射器一致！以下是相关存储器的层级结构：

二、共享存储器

        我们可在共享存储器中声明需要共享的变量，用__shared__限定符，可以和__device__限定符一起使用（可选），一个声明为__shared__的变量有以下特点：

        一、所属的作用域为当前线程块，也即块内的所有线程都都可访问这一限定符所声明的变量，从而起到了块内线程间的数据共享和通信作用；

        二、生命周期和线程保持一致；

        用法：

                 1、__shared__ float array[10];

                       声明了一个大小为10的float数组，块内的线程都共享数据，另外，该数据的大小是固定的，为10

                 2、 extern __shared__ float shared[];

                      声明了一个动态的数据，其大小并没有在此指定，通过执行配置的第三个参数Ns指定（表明共享存储器的空间大小，单位字节）

    

三、利用共享存储器进行矩阵乘法

        除了上面提到的，用__shared__声明的变量可在块内线程间起到通信作用和数据共享外，另外一个作用就是可降低对全局存储器的访问频率，从而节省带宽，当然得存在这样的情况：对全局存储器访问频繁，且能够用共享存储器去代替全局存储器。

        这里以共享存储器为基础，进行了矩阵的乘法运算，这里进行的是矩阵的分块计算，其算法为：对矩阵进行分块为n*n的大小，每个线程块负责计算一个块大小的数据，而具体到每个线程就只计算对应的数据即可。如下图计算C = A *B

            如果不对矩阵进行分块，就是简单的计算，对于计算矩阵C中的每个元素，都要向矩阵A中读取A.width个元素、向矩阵B中读取B.height个元素，如果每个线程负责计算一个对应的数据，由于并发，在计算期间对全局存储器的访问是多么的繁忙，每个线程都要读取A.width + B.height个元素！如果分块计算，那么对于计算一个块内数据，需要向矩阵A读取B.width / BLOCK_SIZE次，向矩阵B读取A.height / BLOCK_SIZE次。通过增加次数来降低对全局存储器的访问频率，节省带宽。

    

四、代码实现

#include <iostream>
#include <cuda_runtime.h>
#include <cstdlib>
#include <cassert>

using namespace std;

//下面定义矩阵分块的信息，用于计算矩阵C，C = A * B
//注意为了处理简单设置A和B的行数和列数相等

const int g_blockSize = 16;		//定义每块的大小
const int g_numBlock = 100;		//块数
const int g_matrixSize = g_blockSize * g_numBlock;	//矩阵的大小

struct Matrix
{
	int width;			//矩阵宽度
	int height;			//矩阵高度
	int stride;			//每行的步长，方便操作子矩阵,非子矩阵：width = stride
	float * pElements;	//数据元素的首地址
};

//创建一个矩阵matrix，其大小为size*size的方阵
//随机数初始化
void InitHostMatrix(Matrix & matrix,int size)
{
	matrix.width = size;
	matrix.height = size;
	matrix.stride = size;
	matrix.pElements = new float[matrix.stride * matrix.height];
	int n = matrix.stride * matrix.height;
	for (int i = 0; i < n; ++i)
	{
		matrix.pElements[i]= 1.f * rand() / RAND_MAX;
	}
}

//销毁一个host矩阵
void DestroyHostMatrix(Matrix & matrix)
{
	delete[] matrix.pElements;
	matrix.pElements = nullptr;
}

//用一个host矩阵创初始化device矩阵
//参数copyContent决定了是否将host矩阵的内容拷贝到device矩阵中,defalt:true
void InitDeviceMatrix(Matrix& deviceMatrix, const Matrix& hostMatrix,bool copyContent = true)
{
	assert(hostMatrix.width > 0 && hostMatrix.height > 0 && hostMatrix.stride && hostMatrix.pElements);
	deviceMatrix.width = hostMatrix.width;
	deviceMatrix.height = hostMatrix.height;
	deviceMatrix.stride = hostMatrix.stride;

	int totalSize = sizeof(float) * hostMatrix.stride * hostMatrix.height;
	cudaError err = cudaMalloc(&deviceMatrix.pElements, totalSize);
	if (err != cudaSuccess)
	{
		std::cerr << "call cudaMalloc fail!" << std::endl;
	}

	if (copyContent)
	{
		err = cudaMemcpy(deviceMatrix.pElements, hostMatrix.pElements, totalSize, cudaMemcpyHostToDevice);
		if (err != cudaSuccess)
		{
			std::cerr << "call cudaMemcpy fail!" << std::endl;
		}
	}
}

//销毁一个Device矩阵
void DestroyDeviceMatrix(Matrix & deviceMatrix)
{
	if (deviceMatrix.pElements)
	{
		cudaFree(deviceMatrix.pElements);
		deviceMatrix.pElements = nullptr;
	}
}

//获取矩阵的子矩阵,根据blockSize来划分块
//声明为__device,表明该函数由设备负责调用
//这里的row、col单位为块，并不是以具体的元素为单位
__device__ Matrix GetSubMatrix(const Matrix matrix, int row, int col)
{
	Matrix resultMatrix;
	resultMatrix.width = g_blockSize;
	resultMatrix.height = g_blockSize;
	resultMatrix.stride = matrix.stride;	//！！！这里设置为父矩阵的步长，以便用GetElement函数能够获取到正确的元素
	resultMatrix.pElements = matrix.pElements + row * g_blockSize * matrix.stride + col * g_blockSize;
	return resultMatrix;
}

//获取矩阵中指定行、列的元素
//这里row、col单位为元素个数，不为块
__device__ float GetElement(const Matrix matrix,int row,int col)
{
	return *(matrix.pElements + matrix.stride * row + col);
}

//设置矩阵指定行、列的元素
__device__ void SetElment(Matrix matrix, int row, int col,float value)
{
	*(matrix.pElements + matrix.stride * row + col) = value;
}


//计算matC,matC = matA * matB
//每个线程负责计算matC中的一个元素，唯一
__global__ void MatMulKernel(const Matrix matA, const Matrix matB, Matrix matC)
{
	//当前计算矩阵块编号
	const int blockRow = blockIdx.y;	//注意x、y是要颠倒的
	const int blockCol = blockIdx.x;

	const int threadRow = threadIdx.y;	
	const int threadCol = threadIdx.x;

	//定义共享子矩阵,每个线程负责从matC中复制数据
	__shared__ float sharedSubA[g_blockSize][g_blockSize];		
	__shared__ float sharedSubB[g_blockSize][g_blockSize];

	Matrix subC = GetSubMatrix(matC, blockRow, blockCol);

	float result = 0.f;
	
	//矩阵matA、matB都被分成了g_numBlock*g_numBlock块
	for (int b = 0; b < g_numBlock; ++b)
	{
		//1、获取当前块的数据
		//根据当前块索引，获取matA、matB对应的子矩阵,并将其内容复制到共享存储器中（每个线程只复制一个数据元素）
		Matrix subA = GetSubMatrix(matA, blockRow, blockCol);
		Matrix subB = GetSubMatrix(matB, blockRow, blockCol);
		//将数据复制到共享存储器中
		sharedSubA[threadRow][threadCol] = GetElement(subA, threadRow, threadCol);
		sharedSubB[threadRow][threadCol] = GetElement(subB, threadRow, threadCol);
		//同步线程
		__syncthreads();
		
		//2、计算
		for (int e = 0; e < g_blockSize; ++e)
		{
			result += sharedSubA[threadRow][e] * sharedSubB[e][threadCol];	//累计结果
		}

		__syncthreads();

	}

	//将结果写入矩阵C对应的子矩阵中
	SetElment(subC, threadRow, threadCol, result);
}


int main()
{
	//定义主机上的矩阵
	Matrix A, B, C;
	InitHostMatrix(A,g_matrixSize);		//初始化host矩阵，行列分别为g_matrixSize
	InitHostMatrix(B,g_matrixSize);
	InitHostMatrix(C,g_matrixSize);

	Matrix dA, dB, dC;
	InitDeviceMatrix(dA, A);
	InitDeviceMatrix(dB, B);
	InitDeviceMatrix(dC, C, false);		//不拷贝C的内容到dC中

	dim3 gridSize(g_numBlock, g_numBlock);		//这里预先指定了g_matrixSize = g_numBlock * g_blockSize
	dim3 blockSize(g_blockSize, g_blockSize);
	MatMulKernel << <gridSize, blockSize>> >(dA, dB, dC);
	//MatMulKernel << <gridSize, blockSize, sizeof(float)*g_blockSize * g_blockSize >> >(dA, dB, dC);

	cudaError err = cudaGetLastError();
	if (err != cudaSuccess)
	{
		std::cerr << "call MatMulKernel!" << std::endl;
	}

	//将计算结构从device拷贝到host
	err = cudaMemcpy(C.pElements, dC.pElements, sizeof(float)*dC.stride * dC.height, cudaMemcpyDeviceToHost);
	if (err != cudaSuccess)
	{
		std::cerr << "call cudaMemcpy fail!" << std::endl;
	}



	//销毁矩阵
	DestroyHostMatrix(A);
	DestroyHostMatrix(B);
	DestroyHostMatrix(C);

	DestroyDeviceMatrix(dA);
	DestroyDeviceMatrix(dB);
	DestroyDeviceMatrix(dC);


	cin.get();
	return 0;
}