GPU并行算法读书笔记-chapter3 CUDA线程模型

目录

• 《GPU并行算法》读书笔记-chapter 3 CUDA线程模型
  • SIMD模型
  • SIMT模型
  • kernel函数
  • CUDA线程结构
    • block
    • grid
    • 线程的全局ID
    • CUDA结构与GPU硬件的映射关系
  • CUDA线程设计
  • 执行配置
    • 什么是执行配置
    • 如何设置执行配置参数
    • 重要知识点总结
  • device函数
  • 自动生成kernel
  • CUDA线程同步

《GPU并行算法》读书笔记-chapter 3 CUDA线程模型

SIMD模型

SIMD（Single Instruction Multiple Data），单指令多数据模型是指计算机执行的指令只有一个，而这条指令处理的数据有多条。

在串行算法中，SIMD模型是用一个计算单元（比如CPU核心），通过切换不同的数据得到最终计算结果。而在并行算法中，SIMD算法在每个计算单元中存入相同的计算指令，各个计算单元独立完成所分配的计算数据的计算。

SIMT模型

CUDA采用SIMD模型，但在SIMD的基础上加以改进，用同一指令序列生成CUDA线程，将多个CUDA线程映射到GPU设备的不同CUDA核，处理不同的数据，这种模型被称为单指令多线程（Single Instruction Multiple Thread）模型。

Q：SIMD和SIMT模型的区别是什么？

A：与SIMD相比，SIMT模型的最大不同在于，CUDA线程是轻量级的线程，线程的切换代价更低，因而程序员可以通过CUDA系统向GPU发射远超CUDA核数量的CUDA线程。这种做法的优势在于：一方面，通过CUDA线程的切换掩盖数据读写延时，另一方面，程序设计中无需考虑并行规模收到GPU硬件的限制，简化并行程序设计。

Q：CUDA实现SIMT模型的机制是什么？

A：CUDA采用两种机制来实现SIMT模型，一种是以kernel函数完成单个线程的指令编码，包括给线程指派要处理的数据，这种机制通常在线程数没有超过CUDA核时使用。二是采用特殊的执行配置调用kernel函数，指定执行kernel函数的线程数量和组织形式。至于线程的生成、切换、回收等工作将根据执行配置给出的信息由CUDA系统完成，无需程序员干预，从而极大地简化了GPU并行计算的实现。

kernel函数

kernel函数是CUDA中最重要的概念之一，是一切GPU端执行代码的七点。在CUDA Fortran中，由attributes(global)定义的kernle函数被CUDA系统复制多份，每份发给一个GPU流处理器形成一个CUDA线程，驱动GPU并行完成计算任务。kernel函数具有如下特点：

• kernel函数只能处理常数和GPU端变量（包括单变量和数组），无论这些变量是参数还是临时变量还是通过使用module得到的共享变量，在kernel中都默认是GPU端变量，如果kernel处理的某个变量是CPU端变量，则产生编译错误。
• kernel函数不能有可选参数，即调用kernel时必须给出每一个参数对应的实参。
• kernel函数中的变量不能有save属性，由于kernel函数是并发执行的，使用具有save属性的变量可能会导致不确定的结果。
• kernel不能被contains包含在其他知己子程序中，也不能通过contains包含任何子程序。
• kernel函数是异步执行的，即调用它的程序在调用后不等待它的返回就继续向下执行，因此，kernel函数不能通过return指定额外的返回点，同样，也不能通过entry提供额外的入口。

kernel函数的执行配置

一般的执行配置格式为：

call kernle<<<grdi,block[,sharedmemsize,strm]>>>(par1,par2)

CUDA线程结构

一般地，程序员在kernel设计时要做的主要是两件事：设计线程的行为；根据线程ID分配这些行为的对象（数据），而获取线程ID就与CUDA的另一个重要概念“CUDA线程结构”有关。

和CPU上共享存储并行的线程不同，GPU可以同时执行成百上千个线程，而通过线程切换GPU中并发执行的线程数更是远超此数，对数量如此庞大的每一个线程单独进行管理代价太大甚至不可实现，为此，CUDA采用了block和grid两级线程层次模型。对线程进行分组、分层管理。

block

根据GPU硬件计算能力的不同，一个block最多可以有1024个线程（具体数值可以查询GPU硬件特性得到），这些线程在逻辑上可以按一维、二维或者三维排列。使用什么维度进行排布由程序员自行决定，但是存在以下三个限制：

• 任意一个维度至少为1
• 如果采用三维排列，z方向最大不能超过64
• block内的总线程数不能超过1024

线程在block中的坐标使用另一个预定义的type（dim3）类型的变量threadIdx来描述，需要注意的是CUDA C中threadIdx的x、y、z 分量从0开始编号，而CUDA Fortran是从1开始编号的。

grid

多个同等大小和形状的线程block组成grid，grid里的block在逻辑上也可以按照一维、二维或者三维排列，但稍有不同的地方是y方向和z方向允许的最大值为65536（2¹⁶)，x方向最大为2³¹,这意味着grid内允许的最大block数量是2⁶³,如果每个block内有1024个线程，那么每次kernel调用最多可以产生2⁷³个线程。一般来讲，这个数字已经可以被简单地理解为“无穷大”，即程序员在求解问题时可以简单地认为GPU拥有“无穷多的运算单元”，而无需考虑GPU设备的实际运算单元数。

线程的全局ID

了解线程的grid和block两级结构，最重要的是在GPU计算中由此获得线程的全局ID。任意线程的全局ID也可以分为x、y、z三个分量，其计算公式为：

idx=(blockidx%x-1)×blockdim%x+threadidx%x

idy=(blockidx%y-1)×blockdim%y+threadidx%y

idz=(blockidx%z-1)×blockdim%z+threadidz%z

线程全局ID是为不同线程分配计算数据的基础，在GPU计算中使用非常频繁，其计算公式需要理解并牢记。

此外，在实际贤臣个调度管理时，CUDA其实还有个更低的层级，即32个线程组成一个线程束（warp），但这个层级是隐含的，程序员无权更改，但在程序优化时仍需要考虑到这一点。

CUDA结构与GPU硬件的映射关系

Q:显然，CUDA线程是从软件的角度去理解GPU的计算概念的，那么，这些CUDA thread，block，grid又是如何与GPU硬件产生联系的呢？

A:GPU的基本计算单眼是流处理器（SP），由kernel函数定义的CUDA线程在SP上执行，每个CUDA线程由一个SP执行，但通过线程切换，每个SP可以执行多个CUDA线程；GPU设备中多个SP组成流多处理器（Streaming Multiprocessor，简称SM），CUDA代码产生的block被发射到SM上执行，每个block由一片SM执行，但同样地，通过block的切换，每片SM可以执行多个block；整个GPU设备由多片SM组成，而单次kernel函数调用产生的grid被CUDA系统映射到整个GPU设备上执行。

CUDA线程设计

我们给出一个例子，来比较CPU和GPU算法的不同从而理解线程设计的原理和原则：

假定有两个同等大小的三维数组，每一个元素都是一个三元矢量，现在要求这两个数组中对应位置三元矢量的点积。

CPU

对于这样的问题，串行的算法很简单，即用三重循环遍历三维数组，计算对应位置的两数组元素（三元矢量）的点积：

do k=1,Nz
	do j=1,Ny
		do i=1,Nx
		c(i,j,k)=dot_product(a(1:3,i,j,k),b(1:3,i,j,k))
		end do
	end do
end do

其中dot_product是Fortran的内置函数，功能是计算两个矢量的点积。

GPU

最直观的并行计算思路是：每一个CUDA线程完成一堆三元矢量的点积，当线程数量足够多时，所有的矢量点积任务都能得到执行。确定算法后，剩下的任务是根据线程ID分配计算任务，这里需要完成的计算任务是三维的，因此假定线程也是按三维排列是最直观的。完整的kernel如下：

module VectorDotProduct
	use cudafor
	contains
	attributes(global) subroutine VDot(a,b,c,Nx,Ny,Nz)
		implicit none
		integer,value :: Nx，Ny，Nz
		real,intent(in),device,dimension(1:3,Nx,Ny,Nz) :: a,b
		real,intent(out),device,dimension(Nx,Ny,Nz) :: c
		integer :: i,j,k
		i=(blockidx%x-1)*blockdim%x+threadidx%x
		j=(blockidx%y-1)*blcokdim%y+threadidx%y
		k=(blockidx%z-1)*blockdim%z+threadidx%z
		if(i>nx .or. j>Ny .or. k>Nz) return
		c(i,j,k)=dot_product(a(1:3,i,j,k),b(1:3,i,j,k))
	end subroutine VDot
end module VectorDotProduct

需要注意的是，kernel函数中加上条件判断是为了防止访存越界。另外，kernel中使用了return语句，但没有指定额外的返回点，因此线程只是提前返回而不是跳转到其他返回点，是符合CUDA Fortran对kernel函数定义的强制约定的，但若后续代码中需要使用线程同步函数，则已经执行return的线程无法达到线程同步点，会导致线程同步失败，此时不能使用return语句。

kernel设计的一般过程

1. 获取线程全局ID
2. 根据线程ID指定线程行为。kernel函数中设定的线程行为类似于CPU串行算法中do循环的循环体，因而将CPU程序移植到GPU上时，转变的只是并行计算思维而不是算法本身。

执行配置

什么是执行配置

CUDA中执行配置的概念是指在调用kernel函数时，与普通的subroutine相比，多了一对三重尖括号“<<< >>>"括起来的一个部分，这个部分我们称之为核函数的执行配置，一般的执行配置的格式为：

call kernel<<<dimgrid,dimblock[,sharedmemszie,stream]>>>(parlist)

fortran中dimgrid，dimblock都是cudafor中预定义的type(dim3)类型的变量，其函数原型为：

type dim3
	integer(kind=4) :: x,y,z
end type

dimgrid和dimblock的定义与赋值方法为：

type(dim3) :: dimgrid,dimblock

dimgrid=dim3(1,1,1)
dimblock=dim3(1,1,1)

核函数不仅在调用方式上和普通的subroutine有一定的区别，另外一个重要的不同是kernel的异步执行：调用kernel只是向GPU发布计算任务并提供诸如线程数量与结构，所需调用参数等，并不等待kernel的执行，即便GPU暂时没有足够的资源启动这个kernel，调用它的CPU代码也在成功发布计算任务之后立即继续执行下一条程序语句。因此，kernel调用结束只能标明“GPU已经成功接受任务”，但不能表示计算任务已经完成，甚至不能表示计算任务已经开始。实际计算任务的完成必须显式执行同步命令cudadevicesynchronize或者其他类似具有阻塞功能的命令才能得到保证。执行此类命令时会等待kernel被执行完毕。

如何设置执行配置参数

在分配dimgrid和dimblock的大小时，必须要保证grid和block大小相乘的结果大于并行的粒度。但为了节约硬件资源，我们也应该使用刚好能够满足此条件的线程分配方法来定义grid和block的大小。比如对于下面的循环：

do i=1,n
	a(i)=b(i)+c(i)
end do

由于是一个一维问题，因此我们只需要定义一维的block和grid即可，假设我们定义的block大小为（32,1,1），则根据前面所说的原则，grid的大小应该为：((n+32-1)/32,1,1)。至于为什么不是n/32，是因为计算机整数除法会直接舍弃小数部分的原因，即对于K=M/N，总有KN≤M，因此，为了确保K和N的乘积不小于M，采用K=(M+N-1)/N的方法。

GPU端线程用什么维度分布更合适并无定规，一般可以根据求解问题方便来选择，对于一维问题，用线程的一维分布（grid，block均为一维）最直观，对于二维问题，block既可以按照二维分布（与物理问题一致），也可以按照一维分布（把物理问题看作多个一维问题），但grid需要使用二维分布，同理，对于三维问题，block既可以按照三维分布，也可以二维分布，还可以一维分布，但grid需要使用三维分布。

重要知识点总结

1. 线程的全局ID

   不管最终的执行配置是怎样的，线程全局ID的三个分量总是可以按照如下的方式求得：

   idx=(blockidx%x-1)*blockdim%x+threadidx%x
   idy=(blockidx%y-1)*blockdim%y+threadidx%y
   idz=(blockidx%z-1)*blockdim%z+threadidx%z
   

2. 执行配置的参数范围限制

   实际的线程分布由执行配置中的grid和block共同决定，在不考虑性能优化的前提下，设计者两个type(dim3)变量的取值只有两个约束条件：

   • 不能超过CUDA规定的线程block和 grid各分量取值的限制。

     对于block，其\(x\)和\(y\)分量的取值范围是1-1024，\(z\)分量的取值范围是1-64，\(x\times y \times z\)的取值范围是1-1024，即：

     \[\left\{\begin{matrix} \begin{align} &1 \le blockdim\%x \le 2^{10}=1024\\ &1 \le blockdim\%y \le 2^{10}=1024\\ &1 \le blockdim\%z \le 2^{6}=64\\ &1 \le blockdim\%x \times blockdim\%y \times blockdim\%z \le 2^{10}=1024 \end{align} \end{matrix}\right. \]

     对于grid，其\(x\)分量的取值范围为\(1\sim 2^{31}\) ，\(y\)分量和\(z\)分量的取值范围均为\(1\sim 2^{16}\),即：

     \[\left\{\begin{matrix} \begin{align} &1\le griddim\%x \le 2^{31}\\ &1\le griddim\%y \le 2^{16}\\ &1\le griddim\%z \le 2^{16}\\ \end{align} \end {matrix}\right. \]

   • 总线程数要确保实际计算任务都能得到执行

     可以将由grid和block给定的线程想象为按照以下规则排列的三位方阵：

     \(x\)方向线程数：\(griddim\%x \times blockdim\%x\)

     \(y\)方向线程数：\(grdidim\%y \times blockdim\%y\)

     \(z\)方向线程数： \(griddim\%z \times blockdim\%z\)

     必须要确保实际要开展的计算任务分解后映射到这样的方阵时，每一维都在方阵内，实际计算任务没有填满这个方阵时允许的（会导致浪费），但不能超过这个方阵（超过部分的计算不会得到线程去执行），即：

     \[\left\{\begin{matrix} \begin{align} &griddim\%x \times blockdim\%x \ge N_{x}\\ &griddim\%y \times blockdim\%y \ge N_{y}\\ &grdidim\%z \times blockdim\%z \ge N_{z}\\ \end{align} \end{matrix}\right. \]

     再结合我们前面说到的一个执行配置的特例，我们可以得到一个通用的设置执行配置的公式：

     \[\begin{align} griddim\%x&=(N_{x}+blockdim\%x-1)/blockdim\%x\\ griddim\%y&=(N_{y}+blockdim\%y-1)/blockdim\%y\\ griddim\%z&=(N_{z}+blockdim\%z-1)/blockdim\%z\\ \end{align} \]

     可以看到的是，grid的维度和实际物理问题的维度是对应的，而blcok的维度和大小都相对自由，但是由于需要覆盖到整个计算任务，因此当block维数较低或者block内线程数较少的时候，可能会导致某个方向的griddim过大，超出其范围。另外还需要注意的是，GPU内部的指令调度最小为一个warp或者半个warp（参见对应GPU型号的架构图），因此在设置block的大小时，x，y，z三个维度的乘积应是单个warp（32）的整数倍，以防止出现浪费。

device函数

与CPU程序一样，kernel函数也可以调用子程序，但只能调用attribute(device)属性的子程序。device子程序的一般格式为：

attributes(device) returntype function funcname(parlist) !function定义格式
attributes(device) subroutine subname(parlis) !subroutine定义格式

device子程序存在以下的限制：

• device子程序只能被global子程序或者其他device子程序调用

• device子程序只能处理常数和GPU段变量，无论这些变量是子程序参数还是临时变量或者通过使用module得到的共享变量，在device中都默认是GPU端变量，如果device子程序中出现CPU端变量将导致编译错误。

• device子程序不能有可选参数

• device子程序中的变量不能具有save属性，由于GPU端代码是并行执行的，使用具有save属性的变量可能会导致不确定的结果。

• device子程序不能递归调用

• device子程序虽然是同步执行的（这一点和kernel不一样），但同样不能通过return指定额外的返回点，也不能通过entry提供额外的入口

• device子程序不能被contains包含在其他主机子程序（subroutine或function中），也不能通过contains包含任何子程序。

• device子程序必须被定义在module中，通过use module的方式进行调用，否则会在编译时报错。

自动生成kernel

自动生成kernel在官方文档中的名称为kernel loop directive，其使用场景主要是能够针对tight nested loop（紧密嵌套的循环）快速生成核函数代码。其一般格式为：

!$cuf kernel do[(n)] <<<grid,block[,sharedmemsize,streamid]>>> ！需要指定sharedmemsize
!$cuf kernel do[(n)] <<<grid,block,stream=streamid>>> !不指定sharedmemsize，仅指定streamid

这条!$cuf kernel do 开头的编译指导语句需要放在do循环的前面，其含义为：将紧邻其后的n重嵌套循环的循环体作为CUDA核函数，按grid，block设定的线程结构来组织CUDA线程并发射到GPU设备中完成计算。do后面紧跟的可选参数n为希望展开的do循环嵌套层数，必须确保紧邻其后的do循环嵌套层数不小于n，如果省略这一参数，则默认只展开紧邻这一编译指导语句的一层do循环。

kernel loop directive example 1：

!cuf kernel do(2) <<<(*,*),(32,4)>>>
do j=1,m
	do i=1,n
	a(i,j)=b(i,j)+c(i,j)
	end do
end do

在这个例子中，编译指导语句定义了一个二维的block，大小为32×4，指明了展开循环层数为2，则循环的主体会按照一下线程结构映射到GPU端：

• threadidx%x为1到32的线程被映射到内层的i循环
• threadidx%y为1到4的线程被映射到外层的j循环

grid的形状，声明为（*,*),是指在编译过程中（m，n为常量）或者运行时（m，n为变量）由CUDA系统计算得出，计算的方式和前文线程设计部分给出的公式一致。当然也可以自行指定。

kernel loop directive example 2：

!$cuf kernel do <<<*,256>>>
do j=1,m
	do i=1,n
	a(i,j)=b(i,j)+c(i,j)
	end do
end do

这里和上一个例子的区别在于没有显式地指明n的大小，而n的默认值为1，因此这里的线程结构会被映射为：

• threadidx%x为1到256的线程映射到外层的j循环上
• 内层i循环在每个线程中串行运行

你也许会想到，考虑到fortran列优先的性质，这里如果让外层串行而内层并行效率会更高，那么我们可以做出如下的修改：

!$cuf kernel do(2) <<<(*,*),(256,1)>>>
do j=1,m
	do i=1,n
	a(i,j)=b(i,j)+c(i,j)
	end do
end do

此时block内的线程总数没有发生变化，维度也没有发生变化，但循环和线程的映射关系却产生了改变：

• 外层j循环为串行执行
• threadidx%x为1到256的线程被映射到内层i循环上

编译指导语句的限制

看似简单清晰的编译指导语句当然也存在着诸多的限制，否则就不必大费周章地去自行设计核函数结构了，其遵循的基本准则为：

• 循环中使用的数组（包括左端项和右端项），必须位于设备端内存。也就是说并不是只需要简单地在CPU代码段前加一句编译指导语句即可，还需要提前做好设备端数组的定义，do循环中数组的替换，后续结果的拷出等等工作
• 循环中使用的单变量（包括左端项和右端项）需要位于设备端内存。如果有单变量位于CPU内存，则CUDA系统会在设备端建立它的拷贝
• 如果设置了do循环的参数n，则紧邻其后的do循环至少要有n重紧密嵌套的循环，可以更多但不能更少。
• 被展开的do循环不能有数据依赖，简单的判断准则是，计算顺序与循环变量无关
• 被展开的n重do循环的循环次数必须在进入循环之前就是确定的，不能根据条件变化或者提前结束。
• 被展开的n重do循环不能有跳转语句goto或提前退出循环的exit
• 被展开的循环体内部可以再包含do循环，但不会被展开为CUDA线程，而是作为CUDA线程的一部分，因而其循环次数无需事先确定，其循环体也可以使用exit
• CUDA Fortran支持的设备端代码可用的子程序（包括内部子程序和自定义的device子程序）均可以正常调用，但不能使用诸如syncthreads，atomic等设备代码专用函数。
• 循环中不能有隐含循环以及Fortran90标准支持的数组整体运算。

CUDA线程同步

GPU计算一般采用并行算法，海量CUDA线程并发执行完成特定计算任务，但在并行计算过程中，有时需要对各CUDA线程的执行进度进行“对齐”，比如，a线程需要用到b线程的计算结果时，必须等b线程已经完成了相关的计算，否则，a线程所用计算结果是不正确的。CUDA系统提供了适应不同需求的线程同步库函数。

同一线程block内的线程同步：syncthread

attributes(device) subroutine syncthread()

该调用只能在kernel函数或者device子程序中使用，同一block内的线程将在调用处对齐，等待所有线程都达到调用处再继续向下执行。若有任何一个线程不能到达调用处（比如，遇到return语句提前返回了，或者syncthread调用位于某个条件选择语句中），则会导致程序执行错误。

当前GPU设备的同步：cudadevicesynchronize

integer function cudadevicesynchronize()

该调用只能在CPU代码中使用，CPU进程/线程在调用处等待当前GPU设备有关操作的完成，如果这些操作是正常完成的，返回值为0，否则返回错误信息ID，GPU数据的赋值、同步拷贝命令等隐含这个调用。