CUDA C编程（五）并行性的表现以及避免分支分化

并行性的表现

  为更好地理解线程束执行的本质，将使用不同的执行配置分析下述的sumMatrixOnGPU2D核函数。使用nvprof配置指标，可以有助于理解为什么有些网格/块的维数组合比其他的组合更好。这些练习会提供网格和块的启发式算法，这是CUDA编程人员必备的技能。

二维矩阵求和的核函数如下所示：

__global__ void sumMatrixOnGPU2D(float *A, float *B, float *C, int NX, int NY)
{
   unsigned int ix = blockIdx.x * blockDim.x + threadIdx.x;
   unsigned int iy = blockIdx.y * blockDim.y + threadIdx.y;
   unsigned int idx = iy * NX + ix;

   if(ix < NX && iy < NY)
   {
      C[idx] = A[idx] + B[idx]
   }
}

用nvprof检测活跃的线程数
  在接下来的部分，将使用生成的sumMatrix对块和网格配置执行试验。首先，需要生成一个参考结果作为性能基准。为此，要先测试一组基础线程块的配置，尤其是大小为（32，32），（32，16），（16，32），和（16，16）的线程块，通过用适当的命令行参数调用sumMatrix测试各种线程块配置。在Tesla M2070上输出以下结果：

  比较这些结果可以看到，最慢的性能是第一个线程块配置（32，32）。最快的线程块配置为（32，16）.这样可以推断出，第二种情况比第一种情况有更多的线程块，因此它的并行性更好。这个理论可以用nvprof和achieved_occupancy指标来验证。一个内核的可实现占用率被定义为：每周期内活跃线程束的平均数量与一个SM支持的线程束最大数量的比值。结果总结如下：（注意，如果系统中有多个GPU,可以用–devices命令行选项指挥nvprof从特定的设备中获取配置信息）

  从结果中可以观察到两件事：1.因为第二种情况中的块数比第一种情况多，所以设备就可以右更多活跃的线程束。其原因可能是第二种情况与第一种情况相比有更高的可实现占用率和更好的性能。2.第四种情况有最高的可实现占用率，但它不是最快的，因此，更高的占用率并不一定意味着有更高的性能。肯定有其他因素限制GPU的性能。

用nvprof检测内存操作
  在sumMatrix内核（C[idx] = A[idx] + B[idx]）中有三个内存操作：两个内存加载和一个内存存储。可以使用nvprof检测这些内存操作的效率。首先，用gld_throughput指标检测内核的内存读取效率，从而得到每个执行配置的差异：

  第四种情况中的加载吞吐量最高，第二种情况中的加载吞吐量大约是第四种情况的一半，但第四种情况却比第二种情况慢。所以，更高的加载吞吐量并不定意味着更高的性能。
  接下来，用gld_efficiency指标检测全局加载效率，即被请求的全局加载吞吐量占所需的全局加载吞吐量的比值。它衡量了应用程序的加载操作利用设备内存带宽的程度。结果总结如下：

  从上述结果可知，最后两种情况下的加载效率是最前面两种情况的一般。这可以解释为什么最后两种情况下更高的加载吞吐量和可实现占用率没有产生较好的性能。尽管在最后两种情况下正在执行的加载数量（即吞吐量）很多，但是那些加载的有效性（即效率）是较低的。
  注意，最后两种情况的共同特征是它们在最内层维数中块的大小是线程束的一半。如前所述，对网格和块启发式算法来说，最内层的维度（block.x）应该总是线程束的大小的倍数。

增大并行性
  一个块的最内层维度(block.x)应该是线程束大小的倍数，这样能极大提高加载效率，因为可以提出如下问题：1.调整block.x会进一步增加加载吞吐量吗？2.有其他办法可以增大并行性吗？带着这些问题，我们可以通过测试sumMatrix使用更大范围的线程配置来回答这些问题：

  从这些结果中可以总结出以下规律:1.最后一次的执行配置块的大小（256，8），这是无效的。一个块中线程总数超过了1024个（这是GPU硬件限制）。2.最好的结果是第四种情况，块大小为（128，2）。3.第一种情况中块的配置为(64，2)，尽管在这种情况下启动的线程块最多，但不是最快的配置。4.因为第二种情况中块的配置为（64，4），与最好的情况有相同数量的线程块，这两种情况应该在设备上显示出相同的并行性。因为这种情况相比（128，2）仍然表现较差，所以可以得出这样的结论：**线程块最内层维度的大小对性能起着关键的作用。**5.在所有其他情况下，线程块的数量都比最好的情况少。因为，增大并行性仍然是性能优化的一个重要因素。
  这样我们可能会想，线程块最少的那些示例应该显示出较低的可实现占用率，线程块最多的那些例子应该显示出较高的可实现占用率。这个理论可以用nvprof检测achieve_occupancy指标来验证一下：

  从上面的结果可以得到，第一种情况(64,2)在所有的例子中可实现占用率最低，但它的线程块是最多的，这种情况在线程块的最大数量上遇到了硬件限制。第四种情况（128，2）和第七种情况(256，2)拥有追随搞得性能配置，有很接近的可实现占用率。在这两种情况下，通过将block.y设置为1来增大块间并行性，观察性能将如何变化。这使得每个线程块的大小减小了，引起了更多的线程块被启动来处理相同数量的数据。这样做会产生如下结果：

  到目前为止，这些配置能产生最佳的性能。特别是(256，1)的块由于（128，1）。可以使用如下的指令查看可实现占用率、加载吞吐量和加载效率：

  值得注意的是，最好的执行配置既不具有最高的可实现占用率，也不具有最高的加载吞吐量。从这些实验中可以推断出，没有一个单独的指标能直接优化性能。我们需要在几个相关的指标间寻找一个恰当的平衡来达到最佳的总体性能。
总结
  在大部分情况下，一个单独的指标不能产生最佳的性能；与总体性能直接相关的指标或事件取决于内核代码的本质；在相关的指标与事件之间寻求一个好的平衡；从不同角度查看内核以寻找相关指标间的平衡；网格/块启发式算法为性能调节提供了一个很好的起点。

避免分支分化

  有时，控制流依赖于线程索引。线程束的条件执行可能引起线程束分化，这会导致内核性能变差。通过重新组织数据的获取模式，可以减少或避免线程束分化。

并行归约问题
  假设要对一个有N个元素的整数数组求和。使用如下的串行代码很容易实现算法：

int sum = 0;
for(int i = 0; i< N; i++)
   sum += array[i];

  如果有大量的数组元素会怎么样呢》如果通过并行计算快速求和呢？鉴于加法的结合律和交换律，数组元素可以以任何顺序求和。所以可以用以下的方法执行并行加法运算：1.将输入向量划分到更小的数据块中；2.用一个线程计算一个数据块的部分和；3.对每个数据块的部分和再求和得出最终结果。
  并行加法的一个常用方法是使用迭代成对实现。一个数据块只包含一对元素，并且一个线程对这两个元素求和产生一个局部结果。然后，这些局部结果在最初的输入向量中就地保存。这些新值被作为下一次迭代求和的输入值。因为输入值的数量在每一次迭代后会减半，当输出向量的长度达到1时，最终的和就已经被计算出来了。
  根据每次迭代后输出元素就地存储的位置，成对的并行求和实现可以被进一步分为以下两种类型：1.相邻配对：元素与它们直接相邻的元素配对；2.交错配对：根据给定的跨度配对元素。具体实现如下图所示：

  尽管以上代码实现的是加法，但任何满足交换律和结合律的运算都可以代替加法。例如，通过调用max代替求和运算，就可以计算输入向量中的最大值。其他有效运算的例子有最小值、平均值和乘积。在向量中执行满足交换律和结合率的运算，被称为归约问题。并行归约问题是这种运算的并行执行。并行规约是一种常见的并行模式，并且是许多并行算法中的一个关键运算。

并行归约中分化
  下图所示的是相邻配对方法的内核实现流程。每个线程将相邻的连个元素相加产生的部分和。在这个内核里，有两个全局内存数组：一个大数组用来存放整个数组，进行归约；另一个小数组用来存放每个线程块的部分和。每个线程块在数组的一部分上独立地执行操作。循环中迭代一次执行一个归约步骤。归约是在就地完成的，这意味着在每一步，全局内存里的值都被部分和代替。_syncthreads语句可以保证，线程块中的任一线程在进入下一次迭代之前，在当前迭代里每个线程的所有部分和都被保存到了全局内存中。进入下一次迭代的所有线程都使用上一步产生的数值。在最后一个循环以后，整个线程块的和被保存进全局内存中。

__global__ void reduceNeighbored(int* g_idata, int* g_odata, unsigned int n)
{
   //set thread ID
   unsigned int tid = threadIdx.x;
   int *idata = g_idata + blockIdx.x * blockDim.x;

   if(idx >= n)
      return;

   for(int stride = 1; stride < blockDim.x; stride *= 2);
   {
      if((tid % (2 * stride)) == 0)
      {
         idata[tid] += idata[tid + stride];
      }
      
      _syncthreads();
   }
   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}

改善并行归约中的分化
  注意以下条件表达式：if((tid % (2 * stride)) == 0),因为这个语句只对偶数ID的线程为true，所以会导致很高的线程树分化。在并行归约的第一次迭代中，只有ID为偶数的线程必须执行这个条件语句的主体，但是所有的线程都必须被调度。在第二次迭代中，只有四分之一的线程是活跃的，但是所有的线程仍然都必须被调度。通过重新组织每个线程的数组索引来强制ID相邻的线程执行求和操作，线程束就能被归约了。下图展示了这种实现，与原来的先比，部分和的存储位置并没有改变，但是工作线程已经更新了。修改后的内核代码如下：

__global__ void reduceNeighbored(int* g_idata, int* g_odata, unsigned int n)
{
   //set thread ID
   unsigned int tid = threadIdx.x;
   unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;
   int *idata = g_idata + blockIdx.x * blockDim.x;

   if(idx >= n)
      return;

   for(int stride = 1; stride < blockDim.x; stride *= 2);
   {
      int index = 2 * stide * tid;
      if((index < blockDim.x)
      {
         idata[tid] += idata[tid + stride];
      }
      
      _syncthreads();
   }
   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}

  对于一个有512个线程的块来说，前8个线程束执行第一轮归约，剩下8个线程束什么也不做。在第二轮里，前4个线程束执行归约，剩下12个线程束什么也不做。因此，这样也彻底不存在分化了。在最后五轮，当每一轮的线程总数小于线程束的大小时，分化就会出现。

交错配对的归约
  与相邻配对方法相比，交错配对方法颠倒了元素的跨度。初始跨度是线程束大小的一半，然后在每次迭代中国减少一半，如下图所示，在每次循环中，每个线程对两个被当前跨度隔开的元素进行求和，以产生一个部分和。与上图相比，交错归约的工作线程并没有变化。但是每个线程在全局内存中的加载/存储位置是不同的。交错归约的内核代码如下：

__global__ void reduceInterleaved(int* g_idata, int* g_odata, unsigned int n)
{
    //set thread ID
   unsigned int tid = threadIdx.x;
   unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;
   int *idata = g_idata + blockIdx.x * blockDim.x;

   if(idx >= n)
      return;

   for(int stride = blockDim.x / 2; stride > 0; stride >> 1);
   {
      if((tid < stride)
      {
         idata[tid] += idata[tid + stride];
      }
      
      _syncthreads();
   }
   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}