从cuda函数库cub看如何高效实现数组加法

从cuda函数库cub看如何高效实现数组加法

---

近期，我们测试了了engine中关于数组求和和zdb中关于数组求和的速度，结果我们发现zdb中关于数组求和的速度快了10%左右，那么cub是如何更高效的利用cuda实现数组求和的呢？我们又可以从中获得什么启发呢？

传统数组求和

如果不考虑并行计算的话，数组求和貌似只有一个方法，从头扫描到尾，但是如果考虑到并行，并且是有用几千核算力的GPU实现并行，可想而知加速比不是一般的大啊。
于是一种显而易见的并行算法就出现了，就是分块。比如我有100,000,000个数据，我把它分成1000个线程去做，每个线程计算1,000,000个数据，最后这1000个线程做加法。
很显然我们就得到了下面的算法：

__global__ void sum_test2(double *src, int64_t len, double *sum){
  int64_t pos_start = blockIdx.x * blockDim.x + threadIdx.x;
  int64_t pos_step = blockDim.x * gridDim.x;

  double thread_sum = 0;

  for(int64_t i=pos_start; i<len; i+=pos_step){
      thread_sum += src[i];
  }

  __shared__ double block_sum;
  block_sum = 0;
  atomicAdd(&block_sum,thread_sum);

  __syncthreads();

  if(threadIdx.x==0){
      atomicAdd(sum,block_sum);
  }
}

cuda实现方案很简单，把这一大堆线程分成若干个block，每个block里面有若干个线程，每个线程负责求出部分和，再把求出同一block里面的所有线程执行结果之和，最后求出所有block的总和。为了避免并发冲突，我们统统采用了原子加操作，我们原先的engine中也是采用的同样的做法，在数据量比较大的情况下，加速比非常可观，可以碾压cpu。

cub优化之后的数组求和

我们按照cub函数库提供的接口实现数组求和，在大数据下面发现效率竟然可以提升10%，调用代码如下：

__global__ void sum_test(double *src, int64_t len, double *sum){
    int64_t pos_start = blockIdx.x * blockDim.x + threadIdx.x;
    int64_t pos_step = blockDim.x * gridDim.x;

    double thread_sum = 0;

    for(int64_t i=pos_start; i<len; i+=pos_step){
        thread_sum += src[i];
    }
    typedef cub::BlockReduce<double,BLOCK_SIZE> BlockReduce;
    __shared__ typename BlockReduce::TempStorage temp_storage;
    double block_sum = BlockReduce(temp_storage).Sum(thread_sum);
    __syncthreads();
    if(threadIdx.x==0){
        atomicAdd(sum,block_sum);
    }
}

我们阅读cub源代码之后发现，他更改了同一block内线程求和这一操作，他并没有使用原子加这一方式，而是采用share memory这一方式，单独拿出一个线程来进行计算，这一个小小的操作速度可以提升很多。
仔细想来，在这一个小小的操作中，因为每个线程做加法的数量相差无几，差不多在同一段时间达到原子加语句，那这样势必会涉及到互斥锁（cuda实现），确实不如直接用一个线程去做来的实在啊！
最后，我们没有调用cub函数库，单独实现了这一操作，速度相差无几：

__global__ void sum_test(double *src, int64_t len, double *sum){
    int64_t pos_start = blockIdx.x * blockDim.x + threadIdx.x;
    int64_t pos_step = blockDim.x * gridDim.x;

    double thread_sum = 0;

    for(int64_t i=pos_start; i<len; i+=pos_step){
        thread_sum += src[i];
    }
    __shared__ double block_sum[128];//已经预定义好为128
    __shared__ double block_sum_1;
    block_sum_1 =0;
    block_sum[threadIdx.x] = thread_sum;
    __syncthreads();
    if(threadIdx.x==0){
        thread_sum_add(block_sum,blockDim.x,&block_sum_1);
        atomicAdd(sum,block_sum_1);
    }
}