warp深度解析
逻辑上，CUDA中所有thread是并行的，但是，从硬件的角度来说，实际上并不是所有的thread能够在同一时刻执行，接下来我们将深入学习和了解有关warp的一些本质。
1. Warps & Thread Blocks
warp是SM的基本执行单元。一个warp包含32个并行thread，这32个thread执行于SIMT模式。也就是说所有thread执行同一条指令，并且每个thread会使用各自的data执行该指令。
block可以是1D、2D或者3D的，但是，从硬件角度看，所有的thread都被组织成一维的，每个thread都有个唯一的ID。每个block的warp数量可以由下面的公式计算获得：
WarpPerBlock = ceil(ThreadPerBlock / warpSize)
一个warp中的线程必然在同一个block中，如果block所含线程数目不是warp大小的整数倍，那么多出的那些thread所在的warp中，会剩余一些inactive的thread，也就是说，即使凑不够warp整数倍的thread，硬件也会为warp凑足，只不过那些thread是inactive状态，需要注意的是，即使这部分thread是inactive的，也会消耗SM资源，这点是编程时应避免的。
2. Warp Divergence(warp分歧)
GPU支持传统的、C-style的显式控制流结构，例如if…else,for,while等等。但和CPU对比来说，GPU没有复杂的分支预测。
这样问题就来了，因为所有同一个warp中的thread必须执行相同的指令，那么如果这些线程在遇到控制流语句时，如果进入不同的分支，那么同一时刻除了正在执行的分支外，其余分支都被阻塞了，十分影响性能。这类问题就是warp divergence。
注意，warp divergence问题只会发生在同一个warp中。 下图展示了warp divergence问题：
为了获得最好的性能，就需要避免同一个warp存在不同的执行路径。避免该问题的方法很多，比如这样一个情形，假设有两个分支，分支的决定条件是thread的唯一ID的奇偶性，kernel函数如下：
__global__ void mathKernel1(float *c)
{
int tid = blockIdx.x * blockDim.x + threadIdx.x;
float a, b;
a = b = 0.0f;
if (tid % 2 == 0)
a = 100.0f;
else
b = 200.0f;
c[tid] = a + b;
}
一种方法是，将条件改为以warp大小为步调，然后取奇偶，代码如下：
__global__ void mathKernel2(void)
{
int tid = blockIdx.x * blockDim.x + threadIdx.x;
float a, b;
a = b = 0.0f;
if ((tid / warpSize) % 2 == 0)
a = 100.0f;
else
b = 200.0f;
c[tid] = a + b;
}
通过测试发现两个kernel函数性能相近，到这里你应该在奇怪为什么二者表现相同呢，实际上是因为当我们的代码很简单，可以被预测时，CUDA的编译器会自动帮助优化我们的代码。（稍微提一下GPU分支预测，这里一个被称为预测变量的东西会被设置成1或者0，所有分支都会执行，但是只有预测变量值为1时，该分支才会得到执行。当条件状态少于某一个阈值时，编译器会将一个分支指令替换为预测指令。）因此，现在回到自动优化问题，一段较长的代码就可能会导致warp divergence问题了。
可以使用下面的命令强制编译器不做优化：
$ nvcc -g -G -arch=sm_20 program.cu -o program
3. Resource Partitioning(资源划分)
一个warp的context包括以下三部分：
1 Program counter
2 Register
3 Shared memory
同一个warp执行context切换是没有消耗的，因为在整个warp的生命期内，SM处理的每个warp的执行context都是“on-chip”的。
每个SM有一个32位register集合放在register file中，还有固定数量的shared memory，这些资源都被thread瓜分了，由于资源是有限的，所以，如果thread数量比较多，那么每个thread占用资源就比较少，反之如果thread数量较少，每个thread占用资源就较多，这需要根据自己的需求作出一个平衡。
资源限制了驻留在SM中blcok的数量，不同的GPU，register和shared memory的数量也不同，就像Fermi和Kepler架构的差别。如果没有足够的资源，kernel的启动就会失败。
当一个block获得到足够的资源时，就成为active block。block中的warp就称为active warp。active warp又可以被分为下面三类：
1 Selected warp
2 Stalled warp
3 Eligible warp
SM中warp调度器每个cycle会挑选active warp送去执行，一个被选中的warp称为Selected warp，没被选中，但是已经做好准备被执行的称为Eligible warp，没准备好要被执行的称为Stalled warp.
warp适合执行需要满足下面两个条件：
1 32个CUDA core有空
2 所有当前指令的参数都准备就绪
例如，Kepler架构GPU任何时刻的active warp数目必须少于或等于64个。selected warp数目必须小于或等于4个（因为scheduler有4个？不确定，至于4个是不是太少则不用担心，kernel启动前，会有一个warmup操作，可以使用cudaFree()来实现）。如果一个warp阻塞了，调度器会挑选一个Eligible warp准备去执行。
CUDA编程中应该重视对计算资源的分配：这些资源限制了active warp的数量。因此，我们必须掌握硬件的一些限制，为了最大化GPU利用率，我们必须最大化active warp的数目。
4. Latency Hiding(延迟隐藏)
指令从开始到结束消耗的clock cycle称为指令的latency。当每个cycle都有eligible warp被调度时，计算资源就会得到充分利用，基于此，我们就可以将每个指令的latency隐藏于issue其它warp的指令的过程中。
和CPU编程相比，latency hiding对GPU非常重要。CPU cores被设计成可以最小化一到两个thread的latency，但是GPU的thread数目可不是一个两个那么简单。
当涉及到指令latency时，指令可以被区分为下面两种：
1 Arithmetic instruction
2 Memory instruction
顾名思义，Arithmetic instruction latency是一个算术操作的始末间隔。另一个则是指load或store的始末间隔。二者的latency大约为：
1 10-20 cycle for arithmetic operations
2 400-800 cycles for global memory accesses
下图是一个简单的执行流程，当warp0阻塞时，执行其他的warp，当warp变为eligible时从新执行。
你可能想要知道怎样评估active warps 的数量来hide latency。Little’s Law可以提供一个合理的估计：
NumberofRequiredWarps = Latency * Throughput
对于Arithmetic operations来说，并行性可以表达为用来hide Arithmetic latency的操作的数目。下表显示了Fermi和Kepler相关数据，这里是以（a + b * c）作为操作的例子。不同的算数指令，throughput（吞吐）也是不同的。
因为memory throughput总是以GB/Sec为单位，我们需要先作相应的转化。可以通过下面的指令来查看device的memory frequency：
$ nvidia-smi -a -q -d CLOCK | fgrep -A 3 "Max Clocks" | fgrep "Memory"
以Fermi为例，其memory frequency可能是1.566GHz，Kepler的是1.6GHz。那么转化过程为：
乘上这个92可以得到上图中的74，这里的数字是针对整个device的，而不是每个SM。
有了这些数据，我们可以做一些计算了，以Fermi为例，假设每个thread的任务是将一个float（4 bytes）类型的数据从global memory移至SM用来计算，你应该需要大约18500个thread，也就是579个warp来隐藏所有的memory latency。
Fermi有16个SM，所以每个SM需要579/16=36个warp来隐藏memory latency。
Occupancy(占用率)
当一个warp阻塞了，SM会执行另一个eligible warp。理想情况是，每时每刻到保证cores被占用。Occupancy就是每个SM的active warp占最大warp数目的比例：
们可以使用cuda库函数的方法来获取warp最大数目：
cudaError_t cudaGetDeviceProperties(struct cudaDeviceProp *prop, int device);
然后用maxThreadsPerMultiProcessor来获取具体数值。
grid和block的配置准则：
- 保证block中thread数目是32的倍数
- 避免block太小：每个blcok最少128或256个thread
- 根据kernel需要的资源调整block
- 保证block的数目远大于SM的数目
- 多做实验来挖掘出最好的配置
Occupancy专注于每个SM中可以并行的thread或者warp的数目。不管怎样，Occupancy不是唯一的性能指标，当Occupancy达到某个值时，再做优化就可能不再有效果了，还有许多其它的指标需要调节
6、Synchronize(同步)
同步是并行编程中的一个普遍问题。在CUDA中，有两种方式实现同步：
1. System-level：等待所有host和device的工作完成
2. Block-level：等待device中block的所有thread执行到某个点
-
cudaDeviceSynchronize
因为CUDA API和host代码是异步的，cudaDeviceSynchronize() 可以用来停下CPU等待CUDA中的操作完成：
cudaError_t cudaDeviceSynchronize(void);
-
synchreads
因为block中的thread执行顺序不定，CUDA提供了一个函数来同步block中的thread。
__synchreads() 函数可以确保同一线程块内的所有线程保持同步，但是不能确保不同线程块直接的线程同步。一个warp内的线程不需要同步；调用一次__synchreads() 至少需要四个时钟周期，一般需要更多时钟周期，应尽量避免使用。
-
memory fence
不保证所有线程运行到同一位置，只保证执行memory fence函数的线程生产的数据能够安全得被其他线程消费；
__threadfence() ：一个线程调用该函数后，该线程在该语句前对全局存储器或者共享存储器的访问已经全部完成，执行结果对grid中的所有线程可见；
__threadfence__block() : 一个线程调用该函数后，该线程在该语句前对全局存储器或共享存储器的访问已经全部完成，执行结果对block中所有线程可见；
7. Bank Conflict
对于同一个wrap中的线程(一个wrap内包含了32个线程)，访问共享存储器时，以half-wrap的形式分两次访问。同一half-wrap内的线程同时可以访问不同的bank，而不同线程对同一个bank 的访问只能顺序进行。
所谓的bank-conflict，就是同一half-wrap内的线程，访问了同一bank里的共享内存。bank-conflict会让原本并行的对共享内存的访存操作变成串行从而极大的降低程序效率。 特殊情况是：half-wrap内所有的线程访问同一个共享内存中的同一地址，会产生一次广播，在这种情况下不会发生bank conflict。
下面有一些小技巧可以避免bank conflict 或者提高global存储器的访问速度
1. 尽量按行操作，需要按列操作时可以先对矩阵进行转置
2. 划分子问题时，使每个block处理的问题宽度恰好为16的整数倍，使得访存可以按照 s_data[tid]=i_data[tid]的形式进行
3. 使用对齐的数据格式，尽量使用nvidia定义的格式如float3,int2等，这些格式本身已经对齐。
4. 当要处理的矩阵宽度不是16的整数倍时，将其补为16的整数倍，或者用malloctopitch而不是malloc。
5. 利用广播，例如s_odata[tid] = tid%16 < 8 ? s_idata[tid] : s_idata[15];会产生8路的块访问冲突而用：s_odata[tid]=s_idata[15];s_odata[tid]= tid%16 < 8 ? s_idata[tid] : s_odata[tid]; 则不会产生块访问冲突
