关于cuda实现归约求和算法的问题

      归约算法的基本思想是，对一个输入数组执行某种计算，然后产生一个更小的结果数组。当大量的数进行加和计算时，可以利用归约算法，多线程同时进行求和计算，使得时间复杂度下降。

算法思想如下：

      设数据总数为N，输入数组为a[N]。首先启用N/2个线程，对于第i个线程，计算 a[i]=a[i]+a[i+N/2]，此时，每个线程都将两个值合并为一个值，得到部分和数组，即a[0]到a[N/2-1]；然后，对部分和数组执行操作，第i个线程计算a[i]=a[i]+a[i+N/4]，部分和数组长度又减半；接着，以此类推，在此操作执行log₂ N次后，a[0]即为数组总和。

       这可通过下面的并行算法实现（记为代码段00）：

输入：a[N]，N=2k                

输出：数组值的总和存储在a[0]中。

BEGIN
  1. p=N/2
  2. While p>0 do
  3. For i=1 to p do in parallel
  4. a[i]=a[i]+a[i+p]
  5. End parallel
  6. p=[p/2]
  7. End while
END

        while循环重复log₂ N次，时间复杂度为O(log₂ N)，与串行时时间复杂度O(N)相比，减小了计算时间。 而此时，Watch out！！

1.不要让“拖后腿的短板”影响程序的正确性

      由于在设备端上各线程的计算速度不同，在执行中，速度快的线程可能会读取速度慢的线程还未写好的那块内存空间。这很可能是很危险的操作。此时，我们可以选择用共享内存，利用同一线程块的线程的同步来抑制这种情况。

      相应的代码段（代码段01）如下：

__shared__ int cache[N];
int tid = threadIdx.x;
cache[tid]=a[tid];//a[]为函数的参数，传入输入数组 __syncthreads();//线程同步
int i=N/2;
while(i!=0)
{
 if (tid<i){
    cache[tid]+=cache[tid+i];
}
__syncthreads();
i=i/2;
}

      Ta-da！记得用共享内存与其同步机制哦，如题所言，不要让拖后腿的慢线程影响你程序的正确性~    

2.__syncthreads()只适用于同一线程块里；not 线程块间

      当数组长度N超过一个线程块内的最大线程数(threadsPerBlock)时，要同时启用多个线程块。相应地将输入数组分成多个组，每组为threadsPerBlock个，每个线程块做上述代码段中的操作，求得部分和（即每组数据的总和）。其中，每个线程块只将threadsPerBlock个（而非N个）数据拷入到共享内存数据，所以数组索引改为以下形式（代码段02）：

if（tid+blockIdx.x*blockDim.x<N）
{
   cache[tid]=a[tid+blockIdx.x*blockDim.x];
}

     求得了部分和（即每组数据的总和），存放在每个线程块的cache[0]里，通过以下语句（代码段03）将每个线程块内计算出的和放入一个数组中：

if（tid==0）
{
partial_a[blockIdx.x]=cache[0];
}

     接着，对此部分和数组也执行代码段01中的操作，进行同样的归约求和，来求得数据总和。

     然而，却不可直接在代码段03下进行数组partial_a[]求和操作，因为每个线程块的执行速度又不同，直接在后面进行归约求和操作，容易发生竞态条件。因为线程块间没有同步机制（这一点值得注意）。此时，可以再调用一个核函数totalAdd()来求总和。

partial_add<<<BlockPerGrid,threadsPerBlock>>>(partial_a); //归约求和里数组和又存放在本数组中,因此两核函数参数可以相同。都为a或partial_a。 

total_add<<<1,BlockPerGrid>>>(partial_a); //取partial_a名字是为形象起见，可把本页partial_a都写为a

3.小心访问到不该访问的内存区域

     以上的讨论限定N=2^k（在代码段00中），若N不是2^k，或者N/threadsPerBlock不是整数，此时需要多开一线程块，而这线程块的线程没有全部用上。cache[threadsPerBlock]后部分有可能没有被赋值为我们的输入数据。则代码段01中两数相加cache[tid]+=cache[tid+i]会使程序变得不正确。给对相应的代码加上限制条件，避免相加这些危险数据。代码如下：

while(i!=0)
{
     if(tid<i)
     {
           if(tid+i+blockIdx.x*blockDim.x<N)
           {
                 cache[tid]+=cache[tid+i];
           }
           else
          {  
                 cache[tid];
          }
      }
      __syncthreads();
      i=i/2;
}

     当然，如果你足够叛逆，不喜欢用前一半数加上后一半数的算法，亦可采用其他算法。比如，相邻两数相加，其和存在前一半数中。更改代码段00：

//a[i]=a[i]+a[i+p]
 a[i]=a[2*i-1]+a[2*i]

      不过，其他地方也需要作出相应修改。

     Ta-da！至此，就可以完成任意大量数据的归约求和算法咯。