本文原载于我的主页：planckscale.info，转载于此。
版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。
上一篇中谈到了编程模型中的Block等概念如何映射到硬件上执行，以及CUDA如何用并行来掩盖延迟。这一篇继续剖析SIMT，谈一谈控制流分叉，指令吞吐和线程间通讯机制。
虽然我们说warp中的线程类似于SIMD,但事实上它是真正的线程。warp中的每一个thread都有自己的指令地址寄存器，允许它们各自执行不同的任务（控制流分叉）。最简单的，比如一个
[php]
if(threadIdx < 10)
{...}
else
{...}
[/php]
语句，将threadIdx=0...31这一个warp划分成两个分支，各自做不同的事情。这个灵活性以性能为代价，当一个warp中控制流出现分叉时，不同分支的线程会被分组相继执行，直到各分支执行完毕后，控制流重新汇聚成一支（上例中即if语句的结束点）。这种情况下执行单元的利用率较低，因为每个分支执行时都需要关闭其他分支的线程，所以这时一些执行单元是用不到的。
为了尽可能高效的计算，需要约束控制流分叉的出现。除了减少流程控制语句外，还需要注意，并不是只要有流程控制语句就一定会带来控制流分叉。关键是，控制流分叉只是针对同一warp中的线程而言，不同warp的线程原本就是串行化执行的，分叉对其无影响。因此，只有流程控制语句的条件在
同一warp内不一致时，才会有控制流分叉。这样,诸如
[php]
if(threadIdx.x / WARPSIZE < n)
{...}
else
{...}
[/php]
这样的语句是不会有分叉的。当然，更宽松的条件如
[php]
if(blockIdx.x < n)
{...}
else
{...}
[/php]
也不会有分叉。依赖于输入数据的条件如
[php]
if(globalArray[threadIdx.x] < n)
{...}
else
{...}
[/php]
则会带来分叉。
CUDA的指令都是针对一个warp中32个线程的并行指令，因而一条指令需要在每个线程中都被执行完才算执行完毕。对于简单的指令如32位浮点数的加、乘，32位整数的加减等，通常都可以由CUDA Core在一个时钟周期内完成，而每个SM中通常都有不少于32个的CUDA Cores，因而对于一个warp中上述类型的简单指令，就是一个CUDA Core处理一个线程，一周期内就可以执行完毕。而对于一些较复杂的指令，执行单元并不能提供这么高的吞吐率，此时一个warp中32个操作需要在多个周期内串行化处理。
我们可以用单位周期内进行的操作数目N除以32来计算指令的吞吐率。以GM204为例，它的SM中有32*4 = 128个CUDA Cores，32个SFU（特殊函数单元），因而在计算32位浮点加法时具有最高吞吐，一个周期内完成128次操作，单位周期内指令吞吐为128/32 = 4；而计算如sin/cos等超越函数时线程不再一一分配到CUDA Cores上，而是要在32个SFU上计算，单位周期内只能完成32次操作，指令吞吐为1条指令每周期.
指令的吞吐率数据可参考CUDA C Programming Guide中 5.4.1. Arithmetic Instructions，该小节以单位时钟周期每SM上能够进行的操作数的形式给出了各指令的吞吐率。
指令吞吐率是我们进行性能优化的有一个重要指标。通常，影响指令吞吐率的因素除了数值计算操作的复杂度、精确度之外，控制流分叉也是一个贡献因子。这里的原因不难理解，控制流分叉时执行单元的利用率下降，使得单位周期内执行的操作数目下降，从而降低了指令吞吐。
到这里，硬件图景下线程的执行就基本说完了，只剩下一个留到最后的话题：线程间交互。通常，不存在任何相互作用的线程，它们之间才能够以任意的顺序执行，像block。但对于warp这样的线程组，是可能与同一block中其他warp通讯或同步的，这时执行顺序就不能任意。所幸即便在block之内，线程间的交互仍然是较弱的，因而底层可以将block划分成warp来分组串行化执行，遇到交互时再另作处理。我们现在来看看这些交互机制。
线程间交互可以细分为通讯和同步两类。通讯主要由公共存储区域交换数据来实现，但也不排除像shuffle这样的特殊方式存在。
从通讯的粒度来看，可以分为warp内部线程间通讯，block内部线程间通讯，block间通讯，更粗的粒度这里不考虑。block之间的通讯则只能基于global memory，block内部的通讯主要基于shared memory/global memory，warp内部线程间除了可以利用上述所有方式，还有一种特殊的shuffle机制.下面我们以通讯的粒度分类陈述各种通讯的实现方式。
block间通讯通常基于两次kernel发射，一次将通讯数据写入global memory，另一次发射读global memory进行后续处理。这种通讯开销较大，主要来自于global memory访存和kernel发射，所以如果有可能，尽量把任务放在一次kernel发射中完成。
或许有人会问，同一个kernel发射中的两个block具有共同的global memory，是不是也可以利用这个特点来构造同一kernel下block间的通讯呢？通常的答案是no，因为block之间执行顺序不定，很难构造有意义的通讯；但如果要较真，答案是yes，我们真的可以构造一些特殊的block间通讯方式。一个例子如下所示，该实例来自于CUDA C Programming Guide B.5. Memory Fence Functions：
[php]
__device__ unsigned int count = 0;
__shared__ bool isLastBlockDone;
__global__ void sum(const float* array, unsigned int N,
volatile float* result)
{
// Each block sums a subset of the input array.
float partialSum = calculatePartialSum(array, N);
if (threadIdx.x == 0) {
// Thread 0 of each block stores the partial sum
// to global memory. The compiler will use
// a store operation that bypasses the L1 cache
// since the "result" variable is declared as
// volatile. This ensures that the threads of
// the last block will read the correct partial
// sums computed by all other blocks.
result[blockIdx.x] = partialSum;
// Thread 0 makes sure that the incrementation
// of the "count" variable is only performed after
// the partial sum has been written to global memory.
__threadfence();
// Thread 0 signals that it is done.
unsigned int value = atomicInc(& count, gridDim.x);
// Thread 0 determines if its block is the last
// block to be done.
isLastBlockDone = (value == (gridDim.x - 1));
}
// Synchronize to make sure that each thread reads
// the correct value of isLastBlockDone.
__syncthreads();
if (isLastBlockDone) {
// The last block sums the partial sums
// stored in result[0 .. gridDim.x-1]
float totalSum = calculateTotalSum(result);
if (threadIdx.x == 0) {
// Thread 0 of last block stores the total sum
// to global memory and resets the count
// varialble, so that the next kernel call
// works properly.
result[0] = totalSum;
count = 0;
}
}
}
[/php]
代码 1. block间通讯实现数组求和
本代码摘录自 CUDA C Programming Guide B.5. Memory Fence Functions
该例实现一个数组的求和，首先各个block计算部分和，然后由最后一个完成部分和计算的block再把所有的部分和加和出最终结果。block间通过一个位于global memory的变量count通讯，它记录了目前已经完成计算的线程数。这样，最后一个完成部分和计算的block就会发现count的数值为最大线
程id，因此可以判定需要由它自己来完成最后从部分和向总和的计算。
不过，为了更好的软件结构，最好还是避免同一kernel的block间产生耦合。同一kernel中block的通讯还涉及到CUDA的weakly-ordered内存模型问题，一个线程中先后两次内存操作在另一个线程看来未必能够保持原有顺序，这产生了相当大的复杂性。我们在下文还会提到这一问题。
block内的线程通讯机制较为丰富，尤其是线程同属一个warp时的shuffle机制。shuffle在Kepler后出现，是一种相当快的线程间通讯方式，它允许同属一个warp的线程间可以互相引用彼此的寄存器，比如下例：
[php]
__global__ void bcast(int arg)
{
int laneId = threadIdx.x & 0x1f;
int value;
if (laneId == 0) // Note unused variable for
value = arg; // all threads except lane 0
value = __shfl(value, 0); // Get "value" from lane 0
if (value != arg)
printf("Thread %d failed.\n", threadIdx.x);
}
[/php]
代码 2. shuffle机制实现一个值向整个warp的广播
本代码摘录自 CUDA C Programming Guide B.14. Warp Shuffle Functions
laneId是warp中线程的一个index，有threadIdx对32取余得到。__shfl(value, 0)语句使得各线程能够访问laneId==0这一线程中value的值。
更常用的通讯机制自然是shared memory和global memory了。其中shared memory更快速，在大多数时候是构建高性能CUDA程序的必由之路。这些常识不再赘述。基于shared/global memory的线程间数据交换，一定要注意线程的同步。block中线程的同步由__syncthreads()实现。线程会等待同block中其他线程都执行到这一点，并且__syncthreads()语句之前的所有shared/global memory操作都尘埃落定，保证block内所有线程在__syncthreads()之后都能看到这些操作的结果。
最后谈一下CUDA采用的weakly-ordered内存模型。它导致一个线程中相继执行的两个存储器操作在另一个线程看来未必是一样的顺序。例如：
[php]
__device__ int X = 1, Y = 2;
//thread 0
__device__ void writeXY()
{
X = 10;
Y = 20;
}
//thread 1
__device__ void readXY()
{
int B = Y;
int A = X;
}
[/php]
代码 3. weakly-ordered内存模型示例
本代码摘录自 CUDA C Programming Guide B.5. Memory Fence Functions
这段代码可能产生A=1，B=20这样的结果。原因是有多种可能的，要么thread 1看到的X、Y的写入顺序被颠倒，要么thread 1中读取顺序被颠倒。这种看似相当毁三观的事情确确实实发生在我们的代码背后。在一个线程里两个相继但无依赖的内存操作，其实际完成的顺序可能是不确定的。在这个线程
看来这并没有导致什么不同，因为两个操作无依赖，并不会破坏因果链；但在另一个线程的眼里，它就暴露出来了。
忍不住插句嘴，这简直就是狭义相对论的世界观在计算机世界的翻版：一个参考系的观察者所看到两个类空间隔事件(可以是相继发生但因距离遥远而无因果关联）在另一个参考系中看来是颠倒的，但有因果关联的两事件在所有观察者看来时序都不会改变。好玩吧？
所以，表面的秩序井然背后有着巨大的复杂性怪兽，为了关牢它的笼子，我们需要约束我们的代码，用合适的机制来实现线程间通讯。要保证另一个线程看起来，两组存储器操作具有我们所希望的顺序，需要用 Memory Fence Function. 这里不再涉及，对更多细节感兴趣的同学，请参考CUDA C Programming Guide B.5. Memory Fence Functions等章节。
（未完待续）