prefix scan异构计算实践

key words: work-efficient, GPU，Parallel

概要

本实验从前缀和问题出发，以串行实现为起点，探讨前缀和问题在串行（单核cpu）执行的不足，从而延伸至异构计算，希望发挥gpu大规模并行的优势，提升计算速度。作为横向对比，还将gpu实现的性能和cpu单核以及c++多线程实现进行对比，探讨为何gpu适合解决该问题。

实验原理和设计

问题定义

前缀和问题可以定义如下：

    y0 = x0
    y1 = x0 + x1
    y2 = x0 + x1+ x2
    ...

前缀和问题可以有如下的输入:一个数据序列（可以理解为数组），一个二元的操作符，一个Identity element。
先讨论二元操作符，对于前缀和问题要在GPU计算实现并行策略，首先要满足一个条件，该OP要满足交换律和结合律那么有哪些操作符符合要求呢？
如max，+， min等
Identity element定义为I op a = a, 即该元素与任意一个数据进行op操作得到该数据本身,在本实验中主要用到加法因此易得0+a = a。

串行实现分析

data: 输入的数组 output: 输出的数组
int now = 0;
for i in [0..size){
  now = now + data[i];
  output[i] = now;
}

假设有一个含n个元素的数组，那么该算法就需要执行 O(n) steps，并且执行O(n)次的加法,这对于数据规模很大的时候的计算是不够理想的，那么如果能发挥GPU的并行计算的优势，能否达到加速呢？

并行算法

Navie

假设我们有一个规模为8的数组

input = [1,2,3,4,5,6,7,8]
output = [1,3,6,10,15,21,28,36]
output = [1, 1+2, 1+2+3, 1+2+3+4,.., 1+..+8]

不难观察，如果我们想要并行的计算出output中的每一个结果，其实对于任意的output[i]都可以看作是一个input[0..i]的
reduction的结果

那么性能如何呢？根据上述的定义，首先已知如果有一个大小为n的数组，我们需要O(logn) steps完成归约操作，但是仍然需要O(n)的work complexity。[2]
因此Naive并行的scan算法需要max(log1,log2,..,log(n-1), log(n)) = O(logn)的step complexity和1+2+3+..+n-1+n = O(n^2)的 work complexity。这样看起来虽然减少了step的规模，但是计算量却比串行的时候还要大，这是无法接受的代价,因此称该算法是work-inefficient的。

Work-efficient

该算法由Blelloch在1990年提出[3],先下一个定论，该算法达到了O(logn)的step complexity 同时拥有O(n)的work complexity,这符合我们的预期，是一个值得实现的并行算法，因此该实验主要比较该算法的性能。
主要由两部分组成，第一部分我们通过logn轮的operation计算出部分的结果，第二部分再将各个部分的结果遍历到数组，从而达到in-place构建数组的结果。

up-sweep

down-sweep

实验内容

串行性能与GPUwork-efficient算法性能对比

实验环境

cpu环境

GPU环境

通过设置多组数据从数据大小4MB到2GB设置多个区间进行测试，以所消耗的时间为指标进行评估

时间对比

加速比

数据解析

通过实验我们可以从加速比第二章图可以得到结论：当数据规模超过8MB时，CUDA计算要快于单核的cpu计算，并且保持小规模的加速比例上升，但加速比例不是特别明显，甚至有些缓慢，这是因为CUDA实现方面不是最优，还有很多优化的空间。另外从第一张图可以看到，随着数据的不断增加，GPU处理的时间也在翻倍，这是因为随着数据规模翻倍，从主机内存拷贝到设备内存的时间也在增加，而拷贝内存的时间又占了一定的比重，overhead因此上升。

总结和讨论

本次实验，通过work-efficient算法，将串行的step complexityO(n)降至O(logn),并且在GPU上实现了该算法，取得了最高1.3倍的加速比。
通过异构计算，如果我们可以把原有的串行任务，划分为多个相互依赖性不强的子任务，并且能够应用良好step complexity以及work complexity的并行算法，那么这样的任务可以分配到GPU上完成，加速计算。
同时因为是异构计算，一旦kernel function返回，cpu又可以执行其他计算任务，但要注意GPU的overhead也是重要的评估指标。若计算任务的overhead远大于实际kernel的计算量那么并不能发挥GPU计算密集型任务的优势。

Future Work

将算法以C++11多线程环境实现并与CUDA计算对比，以及优化现在的算法避免bank confilct[4]等问题,进而取得更大的加速比。

参考文献

[1]:https://en.wikipedia.org/wiki/Prefix_sum
[2]:http://www.enseignement.polytechnique.fr/profs/informatique/Eric.Goubault/Cours09/CUDA/SC07_CUDA_5_Optimization_Harris.pdf
[3]:Blelloch, Guy E. 1990. "Prefix Sums and Their Applications." Technical Report CMU-CS-90-190, School of Computer Science, Carnegie Mellon University.
[4]:https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda

附录

串行实现代码：

void cpu_prefix_sum(int* input, int* output, size_t n){
  output[0] = input[0];
  for(size_t i = 1; i < n; i++){
    output[i] = output[i-1] + input[i];
  }
}

kernel:

__global__ void prefix_sum(int *input, int *output, int *block_sum, size_t data_per_block){
  //copy correspoding data from global memory to shared memory
  //block中threadIdx.x的范围在[0,data_per_block/2 - 1]
  const size_t data_offset = blockIdx.x *data_per_block;
  extern __shared__ int block_data[];
  block_data[2*threadIdx.x] = input[data_offset + 2*threadIdx.x];
  block_data[2*threadIdx.x + 1] = input[data_offset + 2*threadIdx.x + 1];
  //make sure all data is ready
  __syncthreads();
  //up-sweep
  //step = [0~log(data_per_block)-1]
  //offset控制步伐
  int offset = 1;
  for(int d = data_per_block >> 1; d > 0; d >>= 1){
    //在读写之前设置barrier
    __syncthreads();
    if(threadIdx.x < d){
      int a = (2*threadIdx.x+1)*offset - 1;
      int b = (2*threadIdx.x+2)*offset - 1;
      block_data[b] += block_data[a];
    }
    offset <<= 1;
  }
  __syncthreads();
  //进入down-sweep之前需要保存block sum和清0,交给第一个线程去做
  if(0 == threadIdx.x){
    block_sum[blockIdx.x+1] = block_data[data_per_block - 1];
    block_data[data_per_block - 1] = 0;
  }

  //进入down_sweep;
  for(int d = 1; d < data_per_block; d <<= 1){
    offset >>= 1;
    __syncthreads();
    if(threadIdx.x <  d){
      int a = (2*threadIdx.x + 1)*offset - 1;
      int b = (2*threadIdx.x + 2)*offset - 1;
      int temp = block_data[b];
      block_data[b] += block_data[a];
      block_data[a] = temp;
    }
  }
  __syncthreads();

  //copy shared memory to global memory
  output[data_offset + 2*threadIdx.x] = block_data[2*threadIdx.x];
  output[data_offset + 2*threadIdx.x + 1] = block_data[2*threadIdx.x + 1];
}