ECE408@UIUC MP5 CUDA C++实现并行计算的Reduction Tree Addition

ECE408@UIUC CUDA C++实现并行计算的Reduction Tree Addition

如何理解Reduction Tree:

一句话概括：用树形结构实现多线程的并行化从而降低时间复杂度。

Reduction Tree的性质：我们知道求和一个长度为N的数列怎么的也得进行N-1次加法，怎么优化都是进行基于加法的交换/结合定律进行优化。

并行计算设备上实现Reduction Tree的优势：我们可以同时计算很多路，所以我们可以结合足够多的简单的加法计算实现很多个相邻元素的相加结合，归并地进行计算

在实际效果上形成了低时间复杂度的一个树形结构：Reduction Tree，总的时间上从O(N)降低到了O(log(N))。

CUDA C++实现

x
int main(){
  /* Some launch code..... */ 
  //@@ Initialize the grid and block dimensions here
  dim3 DimGrid(numOutputElements, 1, 1);
  dim3 DimBlock(BLOCK_SIZE, 1, 1);
    
  //@@ Launch the GPU Kernel here
  total<<<DimGrid, DimBlock>>>(deviceInput, deviceOutput, numInputElements);
    
  cudaDeviceSynchronize();
  cudaMemcpy(hostOutput, deviceOutput, numOutputElements * sizeof(float), cudaMemcpyDeviceToHost);
    /********************************************************************
   * Reduce output vector on the host
   * NOTE: One could also perform the reduction of the output vector
   * recursively and support any size input. For simplicity, we do not
   * require that for this lab.
   ********************************************************************/
  for (ii = 1; ii < numOutputElements; ii++) {
    hostOutput[0] += hostOutput[ii];
  }
  /* Some ending code..... */ 
}

__global__ void total(float *input, float *output, int len) {
  //@@ Load a segment of the input vector into shared memory
  //@@ Traverse the reduction tree
  //@@ Write the computed sum of the block to the output vector at the
  //@@ correct index
  __shared__ float partialSum[2*BLOCK_SIZE];
  unsigned int t = threadIdx.x;
  unsigned int start = 2*blockIdx.x*blockDim.x;
  unsigned int index1 = start + t;
  unsigned int index2 = index1 + blockDim.x;
  partialSum[t] = (index1 < len) ? input[index1] : 0;
  partialSum[blockDim.x + t] = (index2 < len) ? input[index2] : 0;
  for (unsigned int stride = blockDim.x; stride >= 1; stride /= 2){
    __syncthreads();
    if (t < stride){
      partialSum[t] += partialSum[t+stride];
    }
  }
  if (t == 0){
    output[blockIdx.x] = partialSum[0];
  }
}

Trick 1：为什么grid and block dimensions是这样一个值？

因为对于输出来说每一个block最后都会归并求为一个值，这个值最后的计算由这个block的第一个thread完成，最后输出是每一个block输出一个值，所以为了能输出numOutput这么长的array，那么我们必须有numOutput这么多个block。

Trick 2：为什么Shared Memory是2*BLOCK_SIZE？

因为这个做法的原理是尽可能多的利用一个block里的线程数量，我们假设有512个thread，同时可执行512个计算，那么我们最多能同时计算512个pair，那么能处理的最多的数据就是1024。如果超过的话没法一次算完，不足的话计算资源就会浪费。

Trick 3：为什么实现是partialSum[t] += partialSum[t+stride]; 而非相邻相加？

个人理解是方便地址计算，如果是相邻的话也不是不能算，就是计算的时候没那么方便，因为我们知道/= 2和直接对int移位是一致的，这样有助于减少编译后程序在机器上执行地址计算的复杂程度，一定程度上提高性能。示意图和最终效果如下：