CUDA C编程（六）展开循环以及动态并行

展开循环

  展开循环是一个尝试通过减少分支出现的频率和循环维护指令来优化循环的技术。在循环展开中，循环主体再代码中要多次被编写，而不是只编写一次循环主体再使用另一个循环来反复执行的。任何的封闭循环可将它的迭代次数减少或者完全循环。循环体的复制数量被称为循环展开因子，迭代次数就变为了原始循环迭代次数除以循环展开因子。在顺序数组中，当循环的迭代次数在循环执行之前就已经知道时，循环展开是最有效提升性能的方法。考虑下面的代码：

for(int i = 0; i < 100; i++)
{
   a[i] = b[i] + c[i];
}
如果重复一次循环体，迭代次数能减少到原始循环的一半
for(int i = 0; i < 100; i+=2)
{
   a[i] = b[i] + c[i];
   a[i+1] = b[i+1] + c[i+1];
}

  从高级语言层面上来看，循环展开使性能提高的原因可能不是显而易见的。这种提升来自于编译器执行循环展开时低级指令的改进和优化。例如，在前面循环展开的例子中，条件i<100只检查了50次，而在原来的循环中则检查了100次。另外，因为在每个循环中每个语句的读和写都是独立的，所以CPU可以同时发出内存操作。
  在CUDA中，循环展开的意义非常重大。我们的目标仍然是相同的：通过减少指令消耗和增加更多的独立调度指令来提高性能。因此，更多的并发操作被添加到流水线上，以产生更高的指令的内存带宽。这为线程束调度器提供更多符合条件的线程束，它们可以帮助隐藏指令或内存延迟。

展 开 的 归 约
  如下所示，在reduceInterleaved核函数中每个线程块只处理一部分数据，这些数据可以被认为是一个数据块。如果用一个线程块手动展开两个数据块的处理，会怎么样？reduceUnrolling2是reduceInterleaved核函数的修正版：每个线程块汇总了来自两个数据块的数据。这是一个循环分区的例子，每个线程作用于多个数据块，并处理每个数据块的一个元素：

__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];
}


__global__ void reduceUnrolling2(int* g_idata, int* g_odata, unsigned int n)
{
   unsigned int tid = threadIdx.x;
   unsigned int idx = blockIdx.x * blockDim.x * 2 + threadIdx.x; 

   int *idata = g_idata + blockIdx.x + blockDim.x * 2;

   //汇总数据块
   if(idx + blockDim.x < n)
      g_idata[idx] += g_idata[idx + blockDim.x];
   
   _syncthreads();

   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];
}

  因为现在每个线程块处理两个数据块，我们需要调整内核的执行配置，将网格大小减小至一半：reduceUnrolling2<<<grid.x / 2, block>>>(d_idata, d_odata, size)在编译和运行代码之后，核函数的执行速度比原来快3.42倍。然后接下来一个线程块分别处理4个、8个数据块，相应的速度都得到了提升。在一个线程中有更多的独立内存加载/存储操作会产生更好的性能，因为内存延迟可以更好的倍隐藏起来。可以使用设备内存读取吞吐量指标，以确定这就是性能提高的原因。结果总结为：归约的展开测试用例和设备吞吐量之间是成正比的。

展 开 线 程 的 归 约
  _syncthreads()是用于块内同步的。在归约核函数中，它用来确保在线程进入下一轮之前，每一轮中所有线程已经将局部结果写入全局内存中了。然而，要细想一下只剩下32个或者更少线程(即一个线程束)的情况。因为线程束的执行是SIMT(单指令多线程)的，每条指令之后有隐式的线程束同步过程。因此，规约循环的最后6个迭代可以用下述语句来展开(下面线程束的展开避免了执行循环控制和线程同步逻辑)：

if(tid < 32){
  volatile int *vmem = idata;
  vmem[tid] += vmem[tid + 32];
  vmem[tid] += vmem[tid + 16];
  vmem[tid] += vmem[tid + 8];
  vmem[tid] += vmem[tid + 4];
  vmem[tid] += vmem[tid + 2];
  vmem[tid] += vmem[tid + 1];
  
}

  注意变量vmem是和volatile修饰符一起被声明的，它告诉编译器每次赋值时必须将vmem[tid]的值存回全局内存中。如果省略了volatile修饰符，这段代码将不能正常工作，因为编译器或缓存可能对全局或共享内存优化读写。如果位于全局或共享内存中的变量有volatile修饰符，编译器会假定其值可以被其他线程在任何时间修改或使用。因此，任何参考volatile修饰符的变量强制直接读或写内存，而不是简单的读写缓存或寄存器。基于reduceUnrolling8,线程束的展开可以添加到归约核函数中，如下所示：

__global__ void reduceUnrolling8(int* g_idata, int* g_odata, unsigned int n)
{
   unsigned int tid = threadIdx.x;
   unsigned int idx = blockIdx.x * blockDim.x * 8 + threadIdx.x; 

   int *idata = g_idata + blockIdx.x + blockDim.x * 8;

   //汇总数据块
   if(idx + blockDim.x * 7< n)
   {
      g_idata[idx] += g_idata[idx + blockDim.x];
      g_idata[idx] += g_idata[idx + 2 * blockDim.x];
      g_idata[idx] += g_idata[idx + 3 * blockDim.x];
      g_idata[idx] += g_idata[idx + 4 * blockDim.x];
      g_idata[idx] += g_idata[idx + 5 * blockDim.x];
      g_idata[idx] += g_idata[idx + 6 * blockDim.x];
      g_idata[idx] += g_idata[idx + 7 * blockDim.x];
   }
   _syncthreads();

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

   if(tid < 32)
   {
      volatile int *vmem = idata;
      vmem[tid] += vmem[tid + 32];
      vmem[tid] += vmem[tid + 16];
      vmem[tid] += vmem[tid + 8];
      vmem[tid] += vmem[tid + 4];
      vmem[tid] += vmem[tid + 2];
      vmem[tid] += vmem[tid + 1];
   }

   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}

  这个核函数的执行时间比reduceUnrollWarps8<<<grid.x / 8, block>>>(d_idata, d_odata, size)快1.05倍。结果总结如下：syncthreads能减少新的核函数中的阻塞。

完 全 展 开 时 的 归 约
  如果编译时已知一个循环中的迭代次数，就可以把循环完全展开。因为在Fermi或Kepler架构中，每个块的最大线程数都是1024，并且在这些归约函数中循环迭代次数是基于一个线程块维度的，所以完全展开归约是可能的：

__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 + blockDim.x * 7 < n)
   {
      g_idata[idx] += g_idata[idx + blockDim.x];
      g_idata[idx] += g_idata[idx + 2 * blockDim.x];
      g_idata[idx] += g_idata[idx + 3 * blockDim.x];
      g_idata[idx] += g_idata[idx + 4 * blockDim.x];
      g_idata[idx] += g_idata[idx + 5 * blockDim.x];
      g_idata[idx] += g_idata[idx + 6 * blockDim.x];
      g_idata[idx] += g_idata[idx + 7 * blockDim.x];
   }
    _syncthreads();

   if(blockDim.x >= 1024 && tid < 512)
      idata[tid] += idata[tid + 512];
   _syncthreads();
   
   if(blockDim.x >= 512 && tid < 256)
      idata[tid] += idata[tid + 256];
   _syncthreads();

   if(blockDim.x >= 256 && tid < 128)
      idata[tid] += idata[tid + 128];
   _syncthreads();

   if(blockDim.x >= 128 && tid < 64)
      idata[tid] += idata[tid + 64];
   _syncthreads();

   if(tid < 32)
   {
      volatile int *vmem = idata;
      vmem[tid] += vmem[tid + 32];
      vmem[tid] += vmem[tid + 16];
      vmem[tid] += vmem[tid + 8];
      vmem[tid] += vmem[tid + 4];
      vmem[tid] += vmem[tid + 2];
      vmem[tid] += vmem[tid + 1];
   }
   
   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}

模 板 函 数 的 归 约
  虽然可以手动展开循环，但是使用模板函数有助于进一步减少分支消耗。在设备函数上CUDA支持模板函数。如下所示，可以指定块的大小作为模板函数的参数：相比于reduceCompleteUnrollWarps8,唯一的区别是使用了模板参数替换了块大小。检查块大小的if语句将在编译时倍评估，如果这一条件是false,那么编译时它将会被删除，使得内循环更有效率。例如，在线程块大小为256的情况下调用这个核函数，下述iBlockSize >= 1024 && tid < 512语句将永远是false,编译器会自动从执行内核中移除它。该核函数一定要在switch-case结构中被调用。这允许编译器为特定的线程块大小自动优化代码，但这也意味着它只对在特定块大小下启动reduceCompleteUnroll有效：

switch (blocksize)
{
   case 1024:
       reduceCompleteUnroll<1024><<<grid,x / 8, block>>>(d_idata,d_odata,size);
   case 512:
       reduceCompleteUnroll<512><<<grid,x / 8, block>>>(d_idata,d_odata,size);
   case 256:
       reduceCompleteUnroll<256><<<grid,x / 8, block>>>(d_idata,d_odata,size);
   case 128:
       reduceCompleteUnroll<128><<<grid,x / 8, block>>>(d_idata,d_odata,size);
   case 64:
       reduceCompleteUnroll<64><<<grid,x / 8, block>>>(d_idata,d_odata,size);
}

template <unsigned int iBlockSize>
__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 + blockDim.x * 7 < n)
   {
      g_idata[idx] += g_idata[idx + blockDim.x];
      g_idata[idx] += g_idata[idx + 2 * blockDim.x];
      g_idata[idx] += g_idata[idx + 3 * blockDim.x];
      g_idata[idx] += g_idata[idx + 4 * blockDim.x];
      g_idata[idx] += g_idata[idx + 5 * blockDim.x];
      g_idata[idx] += g_idata[idx + 6 * blockDim.x];
      g_idata[idx] += g_idata[idx + 7 * blockDim.x];
   }
    _syncthreads();

   if(iBlockSize >= 1024 && tid < 512)
      idata[tid] += idata[tid + 512];
   _syncthreads();
   
   if(iBlockSize >= 512 && tid < 256)
      idata[tid] += idata[tid + 256];
   _syncthreads();

   if(iBlockSize >= 256 && tid < 128)
      idata[tid] += idata[tid + 128];
   _syncthreads();

   if(iBlockSize >= 128 && tid < 64)
      idata[tid] += idata[tid + 64];
   _syncthreads();

   if(tid < 32)
   {
      volatile int *vmem = idata;
      vmem[tid] += vmem[tid + 32];
      vmem[tid] += vmem[tid + 16];
      vmem[tid] += vmem[tid + 8];
      vmem[tid] += vmem[tid + 4];
      vmem[tid] += vmem[tid + 2];
      vmem[tid] += vmem[tid + 1];
   }
   
   if(tid == 0)
      g_odata[blockIdx.x] = idata[0];
}

动态并行

  到目前为止，所有核函数都是从主机线程中被调用的。GPU的工作负载完全在CPU的控制下。CUDA的动态并行允许在GPU端直接创建和同步新的GPU内核。在一个核函数中在任意点动态增加GPU应用程序的并行性，是一个令人兴奋的新功能。到目前为止，我们需要把算法设计为单独的、大规模数据并行的内核启动。动态并行提供了一个更有层次结构的方法，在这个方法中，并发性可以在一个GPU内核的多个级别中表现出来。使用动态并行可以让递归算法更加清晰易懂，也更容易理解。有了动态并行，可以推迟到运行时决定需要在GPU上创建多少个块和网格，可以动态地利用GPU硬件调度器和加载平衡器，并进行调整以适应数据驱动或工作负载。在CPU端直接创建工作的能力可以减少在主机和设备之间传输执行控制和数据地需求，因为在设备上执行的线程可以在运行时决定启动配置。

嵌 套 执 行
  在动态并行中，内核执行分为两种类型：父母和孩子。父线程、父线程块或父网格启动一个新的网格，即子网格。子线程块或子网络被父母启动。子网格必须在父线程、父线程块或父网格完成之前完成。只有所有的子网格都完成之后，父母才会完成。下图说明了父网格和子网格的适用范围。主机线程配置和启动父网格，父网格配置和启动子网格。子网格的调动和完成必须进行适当的嵌套，这意味着在线程创建的所有子网格都完成之后，父网格才会完成。如果调用的线程没有显式的同步启动子网格，那么运行时保证父母和孩子之间的隐式同步。在下图中，在父线程中设置了栅栏，从而可以与其子网格显式地同步。

  设备线程中的网格启动，在线程块间是可见的。这意味着，线程可能与由该线程启动的或由相同线程块中其他线程启动的子网格同步。在线程块中，只有当所有线程创建的所有子网格完成之后，线程块的执行才会完成。如果块中所有线程在所有的子网格完成之前退出，那么在这些子网格上隐式同步会被触发。
  当父母启动一个子网格，父线程块与孩子显式同步之后，孩子才能开始执行。父网格和子网格共享相同的全局和常量内存存储，但它们有不同的局部内存和共享内存。有了孩子和父母之间的弱一致性作为保证，父网格和子网格可以对全局内存并发存储。有两个时刻，子网格和它的父线程见到的内存完全相同：子网格开始时和子网格完成时。当父线程优于子网格调用时，所有的全局内存操作要保证对子网格是可见的。当父母在子网格完成时进行同步操作后，子网格所有的内存操作应保证对父母是可见的。
  共享内存和局部内存分别对于线程块或线程来说是私有的，同时，在父母和孩子之间不是可见或一致的。局部内存对线程来说是私有存储，并且对该线程外部不可见。当启动一个子网格时，向局部内存传递一个指针作为参数是无效的。

在 GPU 上 嵌 套Hello World
  为了初步理解动态并行，可以创建一个核函数，使其用动态并行来输出“Hello World”，下图说明了用动态并行完成这个核函数构造的嵌套、递归执行。主机应用程序调用父网格，该父网格在一个线程块中有8个线程。然后，该父网格中的线程0调用一个子网格，该子网格中有一半线程，即4个线程。之后，第一个子网格中的线程0再调用一个新的子网格，这个新的子网格中也只有一半线程，即2两个线程，以此类推，直到最后的嵌套中只剩下一个线程。

  实现这个逻辑的内核代码如下所示。每个线程的核函数执行，会先输出“Hello World”。接着，每个线程检查自己是否该停止。如果在这个嵌套层里线程数大于1，线程0就递归地调用一个带有线程数一般的子网格。

__global__ void nestedHelloWorld(int const iSize,int iDepth)
{
   int tid = threadIdx.x;
   printf("Recursion=%d: Hello World from thread %d block %d\n",iDepth,tid,blockIdx,x);
   
   if(iSize == 1)
      return;
   int nthreads = iSize >> 1;

   if(tid == 0 && nthreads > 0)
   {
      nestedHelloWorld<<<1,nthreads>>>(nthreads,++iDepth);
      printf("------>nested execution depth: %d\n",iDepth);
      
   }
}

  动态并行的限制条件：1.只有在计算能力为3.5或更高的设备上才能被支持；2.通过动态并行调用的内核不能在物理方面独立的设备上启动，然而，在系统中允许查询任一个带CUDA功能的设备性能；3.动态并行的最大嵌套深度为24，但是实际上，在每一个新的级别中大多数内核受限于设备运行时系统需要的内存数量。因为为了对每个嵌套层中的父网格和子网格之间进行同步管理，设备运行时要保留额外的内存。

嵌 套 归 约
  归约可以表示成一个递归函数。在CUDA里使用动态并行，可以确保CUDA里的递归归约核函数的实现像在C语言中一样简单。下面列出了带有动态并行的递归归约的内核代码。原始的网格包含许多线程块，但所有嵌套的子网格中只有一个由其父网格的线程0调用的线程块，核函数的第一步是将全局内存地址g_idata转换为每个线程块的本地地址。接着，如果满足停止条件（这是指如果该条件是嵌套执行树上的叶子），结果就被拷贝回全局内存，并且控制立刻返回给父内核中。如果它不是一片叶子内核，就需要计算本地归约的大小，一半的线程执行就地归约。在就地归约完成后，同步线程块以保证所有部分和的计算。紧接着，线程0产生一个只有一个线程块和一个当前线程块一半线程数量的子网格。在子网格被调用后，所有子网格会设置一个障碍点。因为在每个线程块里，一个线程只产生一个子网格，所以这个障碍点只会同步一个子网格。

__global__ void gpuRecursiveReduce(int *g_idata, int *g_odata, unsigned int isize)
{
   unsigned int tid = threadIdx,x;

   int *idata = g_idata + blockIdx.x * blockDim.x;
   int *odata = &g_odata[blockIdx.x];
   if(isize == 2 && tid == 0)
   {
      g_odata[blockIdx.x] = idata[0] + idata[1];
      return;
   }

   int istride = isize >> 1;
   if(istride > 1 && tid < istride)
   {
      idata[tid] += idata[tid + istride];
   }
   _syncthreads();
   if(tid == 0)
   {
      gpuRecursiveReduce <<<1,istride>>>(idata,odata,istride);

      cudaDeviceSynchronize();
   }
   _syncthreads();
}

  对于一个给定的算法，通过使用不同的动态并行技术，可以有多种可能的实现方式。避免大量嵌套调用有助于减少消耗并提升性能。同步对性能与正确性都至关重要，但减少线程块内部的同步次数可能会使嵌套内核效率更高。因为在每一个嵌套层上设备运行时都要保留额外的内存，所以内核嵌套的最大数量可能是受限制的。这种限制的程度依赖于内核，也可能限制任何使用动态并行应用程序的扩展、性能以及其他的性能。