目录
《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(216),x方向最大为231,这意味着grid内允许的最大block数量是263,如果每个block内有1024个线程,那么每次kernel调用最多可以产生273个线程。一般来讲,这个数字已经可以被简单地理解为“无穷大”,即程序员在求解问题时可以简单地认为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设计的一般过程
- 获取线程全局ID
- 根据线程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需要使用三维分布。
重要知识点总结
-
线程的全局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
-
执行配置的参数范围限制
实际的线程分布由执行配置中的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数据的赋值、同步拷贝命令等隐含这个调用。
标签:kernel,并行算法,读书笔记,do,chapter3,线程,CUDA,GPU,block From: https://www.cnblogs.com/yukina/p/16992715.html