在上一篇博客中，我们提到一个kernel
函数会创建一个网格,每个线程通过 kernel 函数 内置的变量来计算其需要处理的数据索引。在那个 kernel 函数中，我们的任务是计算向量加法，因此我们的kernel 函数也是一维布局的。同 时，我们也提到，线程索引是一个最高支持三维的变量，而网格也是最高支持三维变量。换句话说，网格是线程块的三维数组，线程块是线程的三 维数组。这两个三维变量都可以由dim3这一内置类型表示，对于上一篇博客中的向量和kernel
函数，我们的配置参 数是<<<ceil(n/256.0),256>>>，我们同样可以用下面的代码表示:

1 dim3 dimBlock(128,1,1);
2 dim3 dimGrid(ceil(n/128),1,1);
3 vecAddKernel<<<dimGrid,dimBlock>>>(...);

在上篇博客的配置中，我们只指定了第一维，这样 CUDA 就默认了y和 z 维上都的值都是1，在这里，我们显式的为y 和 z 赋值为1。和dimGrid，dimBlock对应 的 kernel 函数内置变量分别是 gridDim 和 blockDim，分别代表网格的维度和线程块的维度。下图中网格的维度是2维，线程块 的维度是3维，那么启动代码的第一个参数是dim3(2,2,1),第二个参数 是 dim3(2,2,3)。在 kernel 函数中,girdDim.x,gridDim.y,gridDim.z 值分别为2,2,1，blockDim.x,blockDim.y,blockDim.z
分别 为2，2，3。

2.简单矩阵乘法 kernel

矩阵乘法是科学计算中最基础，最重要的运算之一，因此提升矩阵运算的效率很重要。下面我们就实现一个简单的矩 阵乘法 kernel 版。首先，如何把矩阵映射到一个个线程呢？我们很自然的想到，既然矩阵是二维布局的，那么我们就让 线程也按二维布局。如下图所示，网格按二维组织，图中x 轴上一共有5个单位的线程块，y 轴上4个单位，一共20个线程块。每个线程块 大小为16∗16，这样每个矩阵的元素都与一个线程相对应。对于任意一个线程，它所对应的元素索
引就是 Col = blockIdx.x * blockDim.x+threadIdx.x,Row = blockIdx.y * blockDim.y+threadIdx.y。图片中的 线程分为四个区域，用1，2，3，4分别标出。其中，在1的线程块中，每个线程都和一个矩阵元素对应，2中的线程块左 侧的线程和元素对应，由于矩阵列数不一定是线程块 x 轴大小 16 的整数倍，最右侧线程块中最右几列的线程则可能没有元素与之对 应。同样，下方的区域3每个线程块最 下面的线程也没有元素与之对应。在区域4中，下侧和右侧的线程都没有元素与之
对应。

确定了线程和元素的对应关系，下面就在 kernel 函数中进行矩阵乘法,具体代码如下：

 1 /*
 2     简单矩阵乘法
 3     计算 C=A*B
 4     numARows: A的行数
 5     numAColumns:A 的列数
 6     numBRows:B的行数
 7     numBColumns:B的列数
 8     numCRows:C的行数
 9     numCColumns:C的列数
10 */
11 __global__ void matrixMultiply(float *A, float *B, float *C, int numARows,
12                                int numAColumns, int numBRows, int numBColumns,
13                                int numCRows, int numCColumns) {
14 
15 	int Row = blockIdx.y*blockDim.y+threadIdx.y;
16 	int Col = blockIdx.x*blockDim.x+threadIdx.x;
17 	//防止坐标越界
18 	if((Row<numARows)&&(Col<numBColumns)){
19 		float Cvalue = 0.0;
20 		for(int i=0;i<numAColumns;i++){
21 			Cvalue += A[Row*numAColumns+i，]*B[Col+i*numBColumns];
22 		}
23 		C[Row*numCColumns+Col] = Cvalue;
24 	}
25 
26 }

第15行-16行分别计算 C 中的坐标，第 18 行的if判断用来跳过超出矩阵界限的线程(如3右侧部分的线程)。在 for 循环中， 访问A,B,C 中的代码采用的是一维索引，这是因为我们的参数都是指向float的指针， 要把二维的矩阵保存在内存中只能用一维的方式保存，读取当然也只能通过一维的方式。C 语言的二维数组是通过行优先的方式 存储为一维数组的。比如说 A[10][8]中，我们要访问第二行第三个元素，可以通过* (A+1 * 10+3)访问，同理，C的第 Row 行，Col 列可以通过 C[Row*numCColumns+Col]访问到。最后，启动
kernel 的代码如下:

1 #define BLOCK_SIZE 16
2 ...
3  dim3 dimGrid(ceil(numCColumns/float(BLOCK_SIZE)),ceil(numCRows/float(BLOCK_SIZE)),1);
4 
5 dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE,1);
6 
7 matrixMultiply<<<dimGrid,dimBlock>>>(deviceA,deviceB,deviceC,numARows,numAColumns,numBRows,numBColumns,numCRows,numCColumns);

这个简单的矩阵乘法和CPU 串行版差别不大，编码很简单，但是其性能却不足。其原因和 CUDA 的并行执行模型、 CUDA 存储器的特点相关。

3.CUDA 并行执行模型和内存模型

CUDA 是以线程块为单位进行资源分配的，这就导致不同的线程块中的线程无法进行同步，同时，不同的线程块可以乱序执行。每个线程块在执 行时都被调度到不同的多核流处理器上(Streaming Multiprocessor,SM)，一个 GPU 有 若干个 SM，每个 SM 能分配若干个线程块，而每 个线程块能包含的线程数也有上限(现在一般是1024)。这样，能同时执行的线程块数就有上限。线程块分配给 SM 后，将被划分为32个线程为 单位的 Wrap，每个 Wrap 是 threadIdx.x 连续的32个线程。第一个从0-31，第二个从32-63，以此类推。每个
SM 按照 SIMD 的方式执行一 个 warp中的所有线程。在某个时刻，一个 warp 中的所有线程都执行相同的指令。如果有分支指令，如果一个 warp 中的所有线程都满足分支，那么将只执行一个分支的指令，如果部分线程满足一个分支，另一部分满足另一分支 ，分支后的指令将在一个 warp 中执行 两次，一次是满足分支的线程执行，一次是不满足的线程执行。我在上一篇中提到，每个 block 包含的线程数一般是32的倍数正是因 为 warp 的大小是32。

现代计算机体系中大量的设计都是为了减少访问存储器的时间，通过分层的设计在速度和花销上寻求平衡。CUDA 作为一个高性能计算框架，它的设计也体现了这个观点。下图是 CUDA 的内存模型:

图中有四种存储器，他们的访问范围和速度如表所示：

全局存储器是访问速度最慢的存储器，但是它具有容量大的优点，要想减少运行时间，必须减少对全局存储器的访问。寄存器是每个线程私有的，容量很小但是速度最快。常数存储器故名思议里面存储 的数据是只读的，尽管速度也比较慢，但是当 GPU 发现读取常数存储器时，会自动将这些内容存入Cache 中，极大的提高了访问速度。共享存 储器是在全局存储器和寄存器间的折中，它既有较快的速度，也有较大的容量，它之所以被称作共享存储器，是因为它的访问范围是一个线程块，能被一个线程块中的所有线程共享。 回顾我们的简单矩阵乘法 kernel函数，访存和计算主要发生在Cvalue
+= A[Row * numAColumns+i] * B[Col+i * numBColumns] 这行代 码内，在这行代码中，我们访问了两次全局存储器获得了A和B 的一个元素，进行了两次计算(乘法和加法)，浮点运算和访存操作的比值是1：1，而 CPU 的计算 速度比访存速度几个数量级，因此整个运算过程的绝大部分时间都用在了访问内存上。要提高矩阵乘法的运算速度，必须减少访存时间，那么必须 使用更快的存储器。寄存器容量太小几乎无法使用，常数存储器是只读的，那么只能使用共享存储器了，但是共享存储器的容量同样有限，不能存储大量数据。通常，我们把数据集划分为多个块，使每个块都小于共享存储器的容量。

4.更快的矩阵乘法

对于矩阵乘法来说，正好有一种良好的算法可以用，这种算法我们在线性代数里都学过，就是分块矩阵乘法。关于分块矩阵乘法，可以查看维基百科-分块乘法。如下图所示,要计算
C 右下角的分块，我们先把 A 和 B 中的蓝色部分保存到共享存储器中，计算这两个矩阵的乘积，然后再加载 A 和 B 的橙色部分，将矩阵运算的结果和上一次计算的结果相加。一直重复这个过程直到所有的分块计算完毕。分块矩阵乘法的
kernel 函数如下:

 1 #define TILE_SIZE 16 //分块的大小
 2 // Compute C = A * B
 3 __global__ void matrixMultiplyShared(float *A, float *B, float *C, int numARows,int numAColumns, int numBRows,int numBColumns, int numCRows,int numCColumns) {
 4 
 5 	//allocate shared memory;
 6 	__shared__ float destA[TILE_SIZE][TILE_SIZE];
 7 	__shared__ float destB[TILE_SIZE][TILE_SIZE];
 8 
 9 	int tx = threadIdx.x;
10 	int ty = threadIdx.y;
11 
12 	//the x,y index of the thread
13 	int Col = tx + TILE_SIZE*blockIdx.x;
14 	int Row = ty + TILE_SIZE*blockIdx.y;
15 
16 	float CVal = 0.0;
17 	for(int t=0;t<(numAColumns-1)/TILE_SIZE+1;++t){
18 		//将分块中的数据加载到共享内存中
19 		if(Row<numARows && t*TILE_SIZE+tx < numAColumns){
20 			destA[ty][tx] = A[Row*numAColumns+t*TILE_SIZE+tx];
21 		}
22 		else{
23 			destA[ty][tx] = 0.0;
24 		}
25 
26 		if(Col<numBColumns && t*TILE_SIZE+ty<numAColumns)
27 			destB[ty][tx] = B[Col+(ty+t*TILE_SIZE)*numBColumns];
28 		else{
29 			destB[ty][tx] = 0.0;
30 		}
31 		__syncthreads();
32 
33 		for(int j =0;j<TILE_SIZE;j++){
34 			CVal += destA[ty][j]*destB[j][tx];
35 		}
36 		__syncthreads();
37 
38 	}
39 	if(Row<numCRows&&Col<numCColumns)
40 		C[Row*numCColumns+Col] = CVal;
41 }

在第6行和第7行，我们声明了两个共享存储器内的二维数组，__shared__是声明共享存储器的特殊标示符。在13，14行我们得到
了当前线程所对应的 C 矩阵元素的坐标。第16行的 Cval 用来暂存计算结果。在17行起的循环中，我们先是将 A,B 矩阵中 对应元素保存到共享存储器中，如果要加载的索引大于A,B 矩阵的边界，则把这个位置的值设为0，这样不影响计算结果。31 行的__syncthreads()函数用于线程间同步。上面我们说过，共享存储器的访问范围是一个线程块，在分块矩阵乘法中，一个 线程块中的所有线程在一次循环中协同把A,B 矩阵的一个分块的分别一个元素存入共享存储器中。然后再进行计算。要确保所 有的元素都被正确的存入共享存储器，我们就需要让各个线程同步。__syncthreads()的作用就是让一个线程块中的线程在这
里同步，当线程块里的所有线程都执行到这条语句后才继续执行下面的代码，这样就保证了所有的数据都已经被存入共享存储 器中，同样的原因，在第36行我们需要等所有的线程都计算完后再进入下一次循环来读入新的分块。 除了用共享存储器代替全局存储器提高了速度以外，分块矩阵乘法还有一个额外的好处。CUDA 全局存储器采用 DRAM 实现，在现代计算机中 通常采用多提交叉存储的方式来提升多个 DRAM 的访问速度，这种访问方式叫做猝发访问(burst)。在猝发 方式下，给一个地址，DRAM 连续几个单元的数据被相继访问，上面提到过，CUDA
中一个 warp 中的线程按 SIMD 的方式执 行，即执行相同的指令，但是处理的数据是不同的，如果应用程序能利用到连续单元的数据，DRAM 的访问速度将成倍提高。在基本 的矩阵乘法中，每次计算分别访问 A和 B 中的元素，B 的元素下标为Col+i * umBColumns，其 中 Col = blockIdx.x * blockDim.x + threadIdx.x，展开后 B 的元素下标就 是blockIdx.x * blockDim.x + threadIdx.x+i * numBColumns,我们可以看到，同一个线程块中的不同线程它们
的 blockIdx 和 blockDim 都相同，不同的是 threadIdx.x,而且 threadIdx.x 是连续增加的，这就是 说，同一个warp的线程读取的内存单元是连续的，这样B中的元素就可以通过猝发方式访问，
A 中的元素 却不具有这个性质。在分块矩阵乘法中，访问 A 的地址索引是Row * numAColumns+t * TILE_SIZE+tx，其中 tx 就是 threadIdx.x，在 同一个线程块中是线性增加，因此可以合并访问。访问 B 的地址是Col+(ty+t * TILE_SIZE) * numBColumns，其中 Col 是 随 threadIdx.x 线性增加的，也可以合并访问。这样就大大减少了访问全局存储器的时间。