关于 CUDA 中 reduction 运算的优化

最近在Udacity上学习有关基于CUDA的并行运算。在学习当中除了Udacity上面所给的资源以外，我本人还查找了很多其他资料。我计划把我所查找到的这些资料写成技术blog，以供以后查阅。

在Unit 3当中，讲到了三种常见的并行运算：reduce、scan和histogram。简单地讲，reduce运算的输入是一个数据集，和一个运算符（二值，结合性），输出为一个数据，其为数据集中的数据依次与运算符相作用后的结果，如求和、求最大值、求最小值都是reduce. Scan运算的输入是一个数组，一个二值运算符和一个identity element. 这个identity element正如0之于加，1之于乘。常见的scan运算为求累积概率。当然我这种说法很不科学，只要能意会就行了。

最常见的reduce算法就是利用了完全二叉树的结构，分两步完成：

对于这一算法的实现，文献[1]以七种算法的逐级演进说明了reduce的优化过程。下面主要讲讲这个优化思想和过程。

首先一种是二叉树结构的直接实现算法，interleaved addressing:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
__global__ void reduce0(int *g_idata, int *g_odata) {
extern __shared__ int sdata[];

    // each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
sdata[tid] = g_idata[i];
__syncthreads();

    // do reduction in shared mem
     for(unsigned int s=1; s < blockDim.x; s *= 2) {
         if (tid % (2*s) == 0) {
             sdata[tid] += sdata[tid + s];
         }
     __syncthreads();
     }

    // write result for this block to global mem
if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

对于红字那一段，高度分支结构会导致性能下降，因此利用strided index将其修改成非分支结构：

1
2
3
4
5
6
7
for (unsigned int s=1; s < blockDim.x; s *= 2) {
    int index = 2 * s * tid;
    if (index < blockDim.x) {
        sdata[index] += sdata[index + s];
    }
    __syncthreads();
}

由于偶数strided index会产生bank conflicts（以后会讲什么是bank conflicts），我们将其改进为sequential addressing:

1
2
3
4
5
6
for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
    if (tid < s) {
         sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
}

Sequential addressing中，红字凸显的部分说明有一半的线程在第一循环迭代的时候就处于空闲状态。为了改进，我们采用first add during load技术。简单地说，我们采用上述过程一半的block，在刚进入线程的时候就进行一次operation，具体实现如下：

1
2
3
4
5
6
// perform first level of reduction,
// reading from global memory, writing to shared memory
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;
sdata[tid] = g_idata[i] + g_idata[i+blockDim.x];
__syncthreads();

我认为将代码优化到这一地步已经足够了，但是令我惊讶的是，文献[1]进行了如下分析：当前代码的带宽仍没有达到GPU的上限，究其原因，由于对于核心运算的不是存储、运算的辅助性指令造成了代码的瓶颈，代码没有达到最大性能。解决的办法就是拆开循环。

分析指出，当s<=32时，有效线程减少至一束（wrap）。一束线程是SIMD（单指令多数据） synchronous的。这就意味着s<=32时不需要__syncthreads()和if (tid < s)。因此，我们可以拆解最后六个循环：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
for (unsigned int s=blockDim.x/2; s>32; s>>=1)
{
    if (tid < s)
    sdata[tid] += sdata[tid + s];
    __syncthreads();
}
if (tid < 32)
{
    sdata[tid] += sdata[tid + 32];
    sdata[tid] += sdata[tid + 16];
    sdata[tid] += sdata[tid + 8];
    sdata[tid] += sdata[tid + 4];
    sdata[tid] += sdata[tid + 2];
    sdata[tid] += sdata[tid + 1];
}

值得注意的是，这减少了所有线程束的无用功，而不仅仅是最后一个线程束。

进一步地，该文献还指出，在已知迭代数的情况下，可以完全拆解循环。在只讨论二的n次方大小的数组时，我们通过模板，给出blocksize，就可以给出一个通用的reduce算法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
if (blockSize >= 512) {
    if (tid < 256) { sdata[tid] += sdata[tid + 256]; } __syncthreads();
}
if (blockSize >= 256) {
    if (tid < 128) { sdata[tid] += sdata[tid + 128]; } __syncthreads();
}
if (blockSize >= 128) {
    if (tid < 64) { sdata[tid] += sdata[tid + 64]; } __syncthreads();
}
if (tid < 32) {
    if (blockSize >= 64) sdata[tid] += sdata[tid + 32];
    if (blockSize >= 32) sdata[tid] += sdata[tid + 16];
    if (blockSize >= 16) sdata[tid] += sdata[tid + 8];
    if (blockSize >= 8) sdata[tid] += sdata[tid + 4];
    if (blockSize >= 4) sdata[tid] += sdata[tid + 2];
    if (blockSize >= 2) sdata[tid] += sdata[tid + 1];
}

函数声明变为：

1
2
template <unsigned int blockSize>
__global__ void reduce5(int *g_idata, int *g_odata)

最后，文献指出，通过算法级联的思想，可以把上文中first add during load变成multiple add.

1
2
3
4
5
6
7
8
9
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize*2) + threadIdx.x;
unsigned int gridSize = blockSize*2*gridDim.x;
sdata[tid] = 0;
while (i < n) {
    sdata[tid] += g_idata[i] + g_idata[i+blockSize];
    i += gridSize;
}
__syncthreads();

注意这个是更改的开头部分。

参考文献：

[1] Optimizing Parallel Reduction in CUDA - Nvidia


留言

▲