CUDA2.1-原理之索引与warp
本小节来自《大规模并行处理器编程实战》第四节，该书是很好的从内部原理结构上来讲述了CUDA的，对于理解CUDA很有帮助，借以博客的形式去繁取间，肯定会加入自己个人理解，所以有错误之处还望指正。
一、块索引与线程索引
CUDA是细粒度的，数据并行的轻量级线程，在启动一个CUDA的一个Kernel函数的时候，就会创建一个线程网格grid，该网格中的所有线程都是执行该kernel函数的，对于kernel的函数调用形式：kernel<<<dimGrid, dimBlk >>>（argument list）；可以看出，一个函数代表一个网格，里面是由块构成，而每个块中就有许多的轻量级线程了。块是cuda的资源调度单位，块之间交互相较于块内较为困难：
如上图所示，一个网格有许多块，每个块有差不多3维的线程形式，其中线程中的具体的索引就是对应着具体的线程计算，这可以用在矩阵操作中。所以一个kernel对应一个网格，那么什么对应每个块和块中的线程呢？是blockIdx （blockIdx.x ，blockIdx.y，blockIdx.z）和threadIdx（threadIdx.x ，threadIdx.y，threadIdx.z），这两个变量是c语言版本的dim3结构其中有3个无符号整数类型的字段：x、y、z，由CUDA内置的变量，对应的值由运行时系统在运行的时候分配，前者就是用来索引网格中块的位置；后者就是用来索引块中的线程的位置。
一个块中的线程是有上限的（在http://blog.csdn.net/shouhuxianjian/article/details/42427285 的最后有所提及），所以如果只用块中的索引来矩阵计算势必是不够的，这里提下，上面kernel的函数形式中dimGrid, dimBlk 都是三维形式，前者是用来说明这个网格中块在x,y,z 三个维度上的分布；而后者是用来说明每个块中线程在x,y,z 三个维度上的分布。但是根据书中p52页所说，网格是一个二维的块数组，所以网格维度参数中第三个字段z通常是默认为1的，（个人：这里说下，这是这本书说的，那时候还挺早的，现在不一定，因为按照cuda的属性返回值看出来不同的显卡有不同的信息，有兴趣的可以自己测试下）：
例子代码如下：
dim3 dimGrid（32,32,1）; dim3 dimBlk(16,16,16);//这是错的，因为16^3大于512（早期显卡），也大于1024（780ti），所以这里需要注意分配 addCUDA<<<dimGrid, dimBlk>>>(参数列表)；
上面的kernel函数自己写的，这里是为了说明<<<>>>中的两个参数，不同维度上最大值请先调取自己显卡的信息核实再决定，这个函数就是告诉运行时系统启动kernel函数了，通过blockIdx和threadIdx,两个的组合就能生成成千上万，上百万的线程用来计算。
上图就是一个矩阵乘法的例子，行向量乘以列向量作为生成向量的一个元素。如果矩阵很大，那么只靠threadIdx来是无法索引的，可以通过下图这种将整个矩阵划分成不同的子块矩阵来计算，刚好就对应了块索引和块中的元素：
下面给出一个例子，因为暂时看到的资料和教程还没有涉及到2维或者3维的操作，就先按照行主排序或者列主排序来将矩阵划分成1维的向量来计算，不过对于矩阵的元素索引还是遵从2维形式的。
__global__ void MatrixMulKernel(float *Md,float *Pd,int Width){
int Row = blockIdx.y*TILE_WIDTH+threadIdx.y;//计算矩阵的第几行
int Col = blockIdx.x*TILE_WIDTH+threadIdx.x;//计算矩阵的第几列
float Pvalue= 0；//临时变量
for(int k = 0;k<Witdh;++k){
Pvalue += Md[Row*Width+k]*Nd[k*Witdh +Col];//循环计算行向量乘以列向量
Pd[Row*Width +Col] = Pvalue;//将值赋值给结果矩阵
}
}
//执行配置参数
dim3 dimGrid(Width/TILE_WIDTH,Width/TILE_WIDTH);
dim3 dimBlk(TILE_WIDTH,TILE_WIDTH);
//启动核函数，其中WIDTH为矩阵的维度，这里默认方形矩阵，TILE_WIDTH是每个子方形矩阵的维度
MatrixMulKernel<<<dimGrid,dimBlk>>>(Md,Nd,Pd,Width)；
二、同步与透明可扩展性
CUDA中有一个用于一个块中所有线程的栅栏同步函数 __syncthreads()来用于协调这个块中的所有线程，该函数是用来进行块中的线程等待的，即假设一个块中有500个线程，kernel中有工作a，工作b，前250个线程完成了工作a，可是剩下250个线程还没完成，这时候使用了该函数，那么前250个线程就需要等待后续的250个线程完成工作a之后在一起整个500个线程接着完成工作b。该函数具有唯一性，不会被异步，假设在 if -else 语句中:
if（argument） __syncthreads()； else __syncthreads()； end
这种情况下无需担心，因为该函数的意思是要么块中所有的线程都按照其中包含这个函数的路径执行，要么都不执行，不然一半停留在if后面，一半停留在else后面，互相等待，岂不是成了死锁了。
正式因为该函数基于块，所以一个块中相互之间的线程等待时间还能够忍受，如果是基于块的相互等待那么就不能忍受了。正是这种相同块需要栅栏同步，不同块无需栅栏同步，所以不同的显卡上程序才有可伸缩性：
如上图，程序刚好划分成8个块，前者设备每次执行2个块，后者每次可以执行4个块，但是都能很灵活的划分任务进行执行。
三、线程分配
启动kernel函数之后，运行时系统会生成相应的线程网格，然后以块为单位把这些线程分配给这些执行资源。这些执行资源组织成多核流处理器（streaming multiprocessor，SM）。书中这里接着说CUDA运行时系统会自动减少每个SM中分配的块的数量，这个情况是发生在任何一种或者多种资源都不能满足一个SM中所有的块运行的时候。也就是在之前的人为分配好块和线程的参数的时候，在执行的时候运行时系统发现每个sm无法使用满运行的所有的块，那么就自动让那几个块休息不工作（个人：结合CUDA2.2讲述的存储器等的限制，应该是假设一个SM最大支持768个线程，而每个块最大支持256个线程，所以SM算下来最少可以支持3个块，但是如果在kernel中设置好了参数，但是发现由于存储器等限制，那么就减少块数，本来是3块的，现在一个SM只能驻留2个块了）。
这里引用《cuda并行程序设计》6.3章节中一段话：“GPU与CPU不同，GPU不使用寄存器重命名的机制，而是致力于为每个线程都分配真实的寄存器，因此当需要上下文切换时，所需要的操作就是将指向当前寄存器组的选择器或者指针进行更新，指向下一个要执行的线程束的寄存器组，因此是零开销。...如果一个内核函数的每个线程需要的寄存器过多，则每个SM中GPU能够调度的线程块的数量就会受到限制，因此总的可以执行的线程数量也会受到限制。开启的线程数量过少会造成硬件无法被充分利用，性能急剧下降，但是开启过多又意味着资源可能短缺，调度到SM上的线程块数量会减少...由于所使用硬件不同，每个SM可供所有线程使用的寄存器空间大小也不同，分分别有8KB，16KB，32KB以及64KB。牢记，每个线程中每个变量都会占用一个寄存器，举个例子，每个SM拥有32KB的寄存器空间，如果每个线程块有256个线程，则对于每个变量是32位，需要占用4个Bytes的情况下，每个线程可以使用32×1024/4/256=32个寄存器。而每个SM也有所谓最大可使用上限的寄存器数量，如果当前线程块上的寄存器数目是允许的最大值时，每个SM只会处理一个线程块”。
一旦一个块分配给了一个SM，那么该块就会被以大小为32的线程warp所划分，warp的个数的多少（个人：英文版是size，觉得应该说的是分配给每个SM中划分的warp的个数，而每个warp是固定的32个线程）是在具体实现的时候指定的，对于程序员来说这个是透明的，因为无需我们关心。这个参数也可以在不同的显卡的属性信息中得知，在SM中，warp才是线程调度的单位，而不是单个的线程，这个概念不属于CUDA的规范，但是却有助于理解和优化在特定CUDA设备上运行的程序的性能。通常cuda会以半个线程束做一次调度，此时可以将一半的线程束的取数据操作合并成一次连续取数据操作，且指令取一次，并广播给这整个线程束（cuda并行程序设计5.5章节：“由于硬件每次只能为一个线程束获取一条指令...在指令执行层，硬件的调度是基于半个线程束，而不是整个线程束”）。块中相邻的Warp的threadIdx的值是连续的，假设第一个waro包含线程0-31，那么第二个就包含线程32-63，以此类推。
上图就是一个SM，这个sm的上面就是三个不同的块，和每个块划分成不同的warp，每个warp中有32个线程。只要给定块的大小和每个sm中块的数量，就能计算每个sm中驻留的warp的数量了，上图中假设每个块有256个线程，那么每个块就有256/32 = 8个warp了，如果每个sm中只有3个块，那么每个sm就有8×3 = 24个warp，假设在G80类型的显卡中每个sm最多驻留768个线程，每个sm最多驻留768/32 = 24个warp。但是如果说每个SM中只有8个sp（streaming processors）（个人：一个sp执行一个warp还是一个线程？应该是执行一个线程，不过是流水形式的串行执行），那么为什么一个sm中会有这么多的warp。这里就在于CUDA的处理器需要高效的执行长延时操作，比如访问全局存储器（时间比访问寄存器长）。如果warp中线程执行一条指令需要等待前面启动的长延时操作的结果（就是该warp需要从全局存储器中提取数值计算），那么就不选择该warp，而是选择另一个不需要等待结果的驻留的warp（这个warp已经得到了自己需要的结果，所以已经无需等待了，可以直接执行了），当多个warp准备执行的时候，采用优先机制选择一个warp执行，这种机制不产生延时的线程先执行，这就是所谓的延时隐藏（latency hiding）。执行这种处于就绪状态的warp不会带来多余的时间开销，被称之为“零开销线程调度（zero-overhead thread scheduling）”。
所以我们在写网格维度和块的维度的时候，最好考虑到这个，将每个块中的线程的数量设置成32的倍数，而且尽可能的先满足一个块中线程的上限，因为块内线程调度相比较与块间调度来说，时间开销更小，更能加速程序运行，也就是更少的粗粒度，更多的细粒度。（该想法是错的，见下面博客）
顺带看了下其他人博客http://blog.csdn.net/mysniper11/article/details/8269776 的相关理解。其中的：
“每一个块内线程数应该首先是32的倍数，因为这样的话可以适应每一个warp包含32个线程的要求，每一个warp中串行执行，这就要求每一个线程中不可以有过多的循环或者需要的资源过多。但是每一个块中如果线程数过多，可能由于线程中参数过多带来存储器要求过大，从而使SM处理的效率更低。所以，在函数不是很复杂的情况下，可以适当的增加线程数目，线程中不要加入循环。在函数比较复杂的情况下，每一个块中分配32或是64个线程比较合适。每一个SM同时处理一个块，只有在粗粒度层面上以及细粒度层面上均达到平衡，才能使得GPU的利用到达最大。“解答了上面错误想法的地方。
