CUDA C/C++ 从入门到入土 一个入土点：__shared__

在NVIDIA夏令营过程中，最难的一课要数shared memory了，也是本文的重点~

1.什么是shared memory？

唯二在GPU芯片上的存储器就是register和shared memory。其中register就是我们在核函数中直接定义的变量存放空间。 如下图所示：

shared memory 有以下特点：

1. 在GPU芯片上，所以数据传输很快
2. 在一个block内部共享数据
   从这两个特点不难看出，shared memory是一个通信快、可以用于block内部信息交流的利器。也就不难推理出，为什么它如此的重要。

2.shared memory 中一些难点：

从博主亲身经历来看，shared memory的难点根源在于它的第二个特点：在一个block内部共享数据。这会让初学者面临以下一些困难：

1. 在我们刚刚有点习惯并行思维,在kernel函数中做了很多练习，习惯于让下达的指令很好的让各thread很好的并行运转，不互相干扰的时候，shared memroy居然开始让我们的thread之间互相通信，这让我们有点无所适从。
2. 你说通信就通信吧，他只在block内部通信？？？在刚开始的时候，老实说博主对thread，block和grid的认识还比较浅薄，对于只在block内部这件事感到很迷茫。

3.以矩阵乘法–shared memory优化案例为载体破除魔咒

3.1 用shared memory优化矩阵乘法的逻辑

M[m][n] * N[n][k] = c[m][k]

1. 我们知道，矩阵相乘的时候，最终c(i,j)处的值为M矩阵的第i行与N矩阵的第j列逐元素进行乘积并累加后得到。

2. 我们可以分别用一个shared memory对a、b矩阵分块，如下图，M和N矩阵中蓝色部分为申请的两个shared memory：
   shared int L[BLOCK_SIZE][BLOCK_SIZE];
   shared int R[BLOCK_SIZE][BLOCK_SIZE];
   

3. 两个shard memory进行相乘，如下图：
   

4. 可以看到，L矩阵的i行与R矩阵的j列逐元素相乘并求和，得到的值，应当是c(i,j)中的一部分，因为L矩阵的i行与R矩阵的j列逐元素相乘并求和这个行为，本质上是M矩阵的第i行与N矩阵的第j列逐元素进行成绩并累加这个行为的一部分。剩下的一部分怎么补充呢？
   我们只需要让M矩阵的shared memory L矩阵朝右滑动，并让N矩阵的shared memory R矩阵朝下滑动，对应行、列相乘，并反复进行，知道到达矩阵边界处即可。

5. 将前面的每个部分进行求和，即得到了c(i,j)处的值。

3.2 代码解读

有重点在注释，请好好食用~

__global__ void gpu_shared_mult_try2(int *a,int *b, int *c, int m,int n, int k)
{

	//每一个thead找到自己在gird中的坐标
    int ty = blockIdx.y * blockDim.y + threadIdx.y; 
    int tx = blockIdx.x * blockDim.x + threadIdx.x;
    
	//申请shared memory空间，为了方便（后会讲），这里直接申请一个BLOCK的尺寸
    __shared__ int L[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ int R[BLOCK_SIZE][BLOCK_SIZE];
    int sum = 0;
    int index_l;
    int index_r;

	/*
	Q1: 为什么是这个条件？
	if (ty <  ((m / BLOCK_SIZE)+1)*BLOCK_SIZE && tx < ((k / BLOCK_SIZE)+1)*BLOCK_SIZE)
	A1 :  如果条件是 （ty < m && tx<k)，应该非常好理解，即屏蔽掉所有矩阵外的线程，
	只有矩阵内的线程参与运作。在这里条件不一样（可以看出比m和k大一点），是因为我们考虑到到，
	如果我们在滑动shared memory的时候，如果没有刚好到达矩阵边界的情况。在这种情况下，
	最后一个shared memory会超过一点点矩阵边缘，为了让超出的部分不影响我们矩阵运算，
	我们需要对超出的部分置0，但是，我们不能直接屏蔽掉超出的部分！！！ 
	置0也需要我们下达操作的命令！ 不能直接屏蔽！因为这部分超出矩阵的线程，
	如果被屏蔽了，他们发生什么事情我们是不可控的！*/
    if (ty <  ((m / BLOCK_SIZE)+1)*BLOCK_SIZE && tx < ((k / BLOCK_SIZE)+1)*BLOCK_SIZE)
    {
    	//我们在每次滑动后，进行矩阵相乘，得到理论中说的每一“部分”
        for (int i=0; i < n/BLOCK_SIZE; i++)
        {	
        	//我们需要首先赋值给shared memory，然后两个shared memory进行相乘
        	//计算坐标应该是基本功啦~，记得画图配合理解。
            index_l = ty * n + i * BLOCK_SIZE + threadIdx.x;
            index_r = (i * BLOCK_SIZE + threadIdx.y) * k + tx;
            /*
            Q2：怎么理解直接赋值给L[threadIdx.y][threadIdx.x] ？
            A2：涉及到我们对shared memory “block内共享”的理解。这里其实有一个巧合，
            对每个thread来说，threadIdx.x和threadIdx.y
            这两个数字都是只在block内部独一无二的，
            一旦跨越到另一个block，threadIdx又会从0开始计数（只不过blockIdx变化了）。
            而shared memory刚好是在block内部共享，
            刚好我们又为它申请到了[BLOCK_SIZE][BLOCK_SIZE]的空间，
            也就是说一整个BLOCK大小，刚好可以装下一个BLOCK内所有的thread。
            所以，我们可以直接在核函数中让每一个线程在BLOCK当中的坐标，
            和shared memory中在BLOCK当中的坐标对应上就好了，
            因此，可以直接赋值给L[threadIdx.y][threadIdx.x] 。
            我们可以形象得理解为，恰好有一个矩阵的钉子挂在天上和一个同样大小排布的矩阵的洞，
            钉子从天上一次落下落到洞里就行，这也是为什么后续需要_syncthreads();进行同步，
            因为你得确保所有钉子都掉洞里了。
            另外，我们在这里还实现了赋0的操作。
            */
            L[threadIdx.y][threadIdx.x] = (i * BLOCK_SIZE + threadIdx.x) < n ? a[index_l]:0;
            R[threadIdx.y][threadIdx.x] = (threadIdx.y + i * BLOCK_SIZE) < n ? b[index_r]:0;
            __syncthreads();
            //随后，把之前的几个“部分”累加即可
            for (int i=0; i < BLOCK_SIZE;i++)
            {
                sum += L[threadIdx.y][i] * R[i][threadIdx.x];
            
            }
            __syncthreads();
        
        }
        //赋值给c
        c[ty * k + tx] = sum;
    }



}