面向GPU计算平台的归约算法的性能优化研究

1 GPU归约算法的实现与优化

图3-1为本文提出的GPU归约算法总图，GPU归约求和算法的实现可以定义为三个层次：

1. 线程内归约：线程从global memory中读取一个或多个数据进行归约操作，再把归约结果写入至LDS；
2. work-group内归约：work-group对LDS的数据进行内部归约操作，求出局部归约结果；
3. work-group间归约：对每一个work-group所得的局部归约结果进行累加操作，得到最终归约结果。

本节将会以Naïve Reduction为起点，逐步地探求并行归约算法的优化要素，以最大化地提升算法性能。

GPU归约算法的Naïve实现采用分治思想，将原始数据划分为多个块；然后对每个块进行局部归约操作，求出块内的局部归约结果，最后再对局部归约结果进行全局归约操作，得到最终归约结果。本文的归约算法优化均以Naïve Reduction为基础进行的。

1.1GPU归约算法的Naïve实现

reduce baseline算法介绍

Baseline算法比较简单，分为三个步骤。第一个步骤是将数据load至shared memory中，第二个步骤是在shared memory中对数据进行reduce操作，第三个步骤是将最后的结果写回global memory中。代码如下：

__global__ void reduce0(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    //each thread loads one element from global memory to shared mem
    unsigned int i=blockIdx.x*blockDim.x+threadIdx.x;
    unsigned int tid=threadIdx.x;
    sdata[tid]=d_in[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)d_out[blockIdx.x]=sdata[tid];
}

GPU归约算法的Naïve实现采用分治思想，将原始数据划分为多个块；然后对每个块进行局部归约操作，求出块内的局部归约结果，最后再对局部归约结果进行全局归约操作，得到最终归约结果。本文的归约算法优化均以Naïve Reduction为基础进行的，其算法伪代码如下：

#include <cuda_runtime.h>
#include <iostream>

__global__ void reduceSumKernel(float *src, float *dest, int n) {
    extern __shared__ float lSum[];
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int idx_loc = threadIdx.x;

    // 线程内归约
    lSum[idx_loc] = (idx < n) ? src[idx] : 0;
    __syncthreads();

    // 线程块内归约
    for (int i = 1; i < blockDim.x; i = 2 * i) {
        if (idx_loc % (2 * i) == 0) {
            lSum[idx_loc] += lSum[idx_loc + i];
        }
        __syncthreads();
    }

    // 线程块间归约
    if (idx_loc == 0) {
        atomicAdd(dest, lSum[0]);
    }
}
int main()
{
    const int N = 1024 * 1024; // 数据大小
    const int blockSize = 256; // 线程块大小
    const int numBlocks = (N + blockSize - 1) / blockSize; // 线程块数量

    // 主机端数据
    float *src;
    float *dest;
    src = new float[N];
    dest = new float[1];

    // 初始化数据
    for (int i = 0; i < N; i++)
    {
        src[i] = 1.0f;
    }
    dest[0] = 0.0f;

    // 设备端内存分配
    float *d_src;
    float *d_dest;
    cudaMalloc(&d_src, N * sizeof(float));
    cudaMalloc(&d_dest, sizeof(float));

    // 数据传输到设备
    cudaMemcpy(d_src, src, N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_dest, dest, sizeof(float), cudaMemcpyHostToDevice);

    // 调用内核
    reduceSumKernel<<<numBlocks, blockSize, blockSize * sizeof(float)>>>(d_src, d_dest, N);
    cudaDeviceSynchronize();

    // 数据从设备传输回主机
    cudaMemcpy(dest, d_dest, sizeof(float), cudaMemcpyDeviceToHost);

    // 输出结果
    std::cout << "Sum: " << dest[0] << std::endl;

    // 验证结果
    float expectedSum = static_cast<float>(N);
    if (dest[0] == expectedSum)
    {
        std::cout << "Result is correct." << std::endl;
    }
    else
    {
        std::cout << "Result is incorrect." << std::endl;
    }

    // 释放内存
    delete[] src;
    delete[] dest;
    cudaFree(d_src);
    cudaFree(d_dest);

    return 0;
}

1.2 GPU归约算法的优化

1.2.1线程内归约优化

线程内归约是归约算法在GPU的移植与优化中常常得不到重视的内容。绝大多数的归约算法的GPU实现和优化都把work-group内归约优化作为算法优化核心，然而，线程内归约才是GPU归约算法影响性能的关键因素，本节对线程内归约过程展开详细的讨论与分析。

Naïve Reduction没有进行线程内归约，一个线程仅仅对应一个数据，仅负责将数据从global memory加载至LDS中，然后在LDS中进行work-group内归约。由于没有进行线程内归约优化，在随之进行的work-group内归约从第一层归约开始，便有一半线程是处于空闲状态，极大地造成了计算资源的浪费。

为了更充分地利用计算资源，应尽可能的使所有线程均参与归约操作，将空闲线程出现的时间尽可能地往后“推移”。因此，在work-group内归约开始之前进行线程内归约操作：每个线程对应多个数据，线程从global memory依次读取多个数据并对其进行归约操作，然后再把归约结果写入LDS。线程内归约将每个线程简单的数据加载操作转变为加载归约操作(把原本每次只加载一个数据变成加载多个数据并归约累加，把累加结果写入LDS中)。这里需要注意的是，我们将每个线程进行线程内归约时处理数据的数目定义为线程内归约粒度。因此，在进行work-group内归约之前，所有线程均参与了归约操作，提升了线程计算量和资源利用率，从而挖掘出归约算法更多的并行潜力。

Global-Stride Kernel

每一个线程以全局线程总数(global stride)为步长，依次读取相距global stride 的多个数据(数据量由线程内归约粒度times控制)，然后对这些数据进行归约处理，最后把归约结果写入到位于LDS中的lSum数组，再进行下一层次的work-group内归约优化。其伪代码如下所示：

#include <cuda_runtime.h>
#include <iostream>

__global__ void reduceSumKernel(float *src, float *dest, int n) {
    extern __shared__ float lSum[];

    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int idx_loc = threadIdx.x;
    int globalSize = gridDim.x * blockDim.x;

    // 线程内归约
    float temp = 0;
    for (int i = 0; i < 1024 && idx + i * globalSize < n; i++) {
        temp += src[idx + i * globalSize];
    }
    lSum[idx_loc] = temp;
    __syncthreads();

    // 线程块内归约
    for (int i = 512; i > 0; i /= 2) {
        if (idx_loc < i) {
            lSum[idx_loc] += lSum[idx_loc + i];
        }
        __syncthreads();
    }

    // 线程块间归约
    if (idx_loc == 0) {
        atomicAdd(dest, lSum[0]);
    }
}

int main() {
    const int N = 1024 * 1024; // 数据大小
    const int blockSize = 1024; // 线程块大小
    const int numBlocks = (N + blockSize - 1) / blockSize; // 线程块数量
    // 主机端数据
    float *src;
    float *dest;
    src = new float[N];
    dest = new float[1];

    // 初始化数据
    for (int i = 0; i < N; i++)
    {
        src[i] = 1.0f;
    }
    dest[0] = 0.0f;

    // 设备端内存分配
    float *d_src;
    float *d_dest;
    cudaMalloc(&d_src, N * sizeof(float));
    cudaMalloc(&d_dest, sizeof(float));

    // 数据传输到设备
    cudaMemcpy(d_src, src, N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_dest, dest, sizeof(float), cudaMemcpyHostToDevice);

    // 调用核函数
    reduceSumKernel<<<numBlocks, blockSize, blockSize * sizeof(float)>>>(d_src, d_dest, N);
    cudaDeviceSynchronize();
   
    // 数据从设备传输回主机
    cudaMemcpy(dest, d_dest, sizeof(float), cudaMemcpyDeviceToHost);
    // 输出结果
    std::cout << "Sum: " << dest[0] << std::endl;

    // 验证结果
    float expectedSum = static_cast<float>(N);
    if (dest[0] == expectedSum)
    {
        std::cout << "Result is correct." << std::endl;
    }
    else
    {
        std::cout << "Result is incorrect." << std::endl;
    }
    // 释放内存
    delete[] src;
    delete[] dest;
    cudaFree(d_src);
    cudaFree(d_dest);
    return 0;
}

Local-Stride Kernel

 每一个线程以work-group内线程数(local stride)为步长，读取相距local stride的多个数据(数据量由线程内归约粒度times控制)，然后对这些数据进行归约处理，最后把归约结果写入到位于LDS中的lSum数组，再进行下一层次的work-group内归约优化。其伪代码如下所示：

 两个kernel的区别在于每个线程读取相邻数据的步长不同。

#include <iostream>
#include <cuda_runtime.h>

__global__ void reduceSumKernel(float *src, float *dest, int data_len) {
    extern __shared__ float lSum[];

    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int idx_loc = threadIdx.x;
    int globalSize = gridDim.x * blockDim.x;

    // 线程内归约
    float temp = 0;
    for (int i = 0; i < 1024 && idx + i * globalSize < data_len; i++) {
        temp += src[idx + i * globalSize];
    }
    lSum[idx_loc] = temp;
    __syncthreads();

    // 线程块内归约
    for (int i = 512; i > 0; i /= 2) {
        __syncthreads();
        if (idx_loc < i) {
            lSum[idx_loc] += lSum[idx_loc + i];
        }
    }

    // 线程块间归约
    if (idx_loc == 0) {
        atomicAdd(dest, lSum[0]);
    }
}

int main() {
    const int N = 1024 * 1024; // 数据长度
    const int blockSize = 1024; // 每个block的线程数
    const int numBlocks = N / blockSize; // block的数量

    float *src, *dest;
    float *d_src, *d_dest;

    // 主机内存分配
    src = new float[N];
    dest = new float[1];

    // 初始化数据
    for (int i = 0; i < N; i++) {
        src[i] = 1.0f;
    }

    // 设定dest为0
    dest[0] = 0.0f;

    // 设备内存分配
    cudaMalloc(&d_src, N * sizeof(float));
    cudaMalloc(&d_dest, sizeof(float));

    // 数据复制到设备
    cudaMemcpy(d_src, src, N * sizeof(float), cudaMemcpyHostToDevice);

    // 调用内核函数
    reduceSumKernel<<<numBlocks, blockSize, blockSize * sizeof(float)>>>(d_src, d_dest, N);

    // 数据复制回主机
    cudaMemcpy(dest, d_dest, sizeof(float), cudaMemcpyDeviceToHost);

    // 输出结果
    std::cout << "Sum: " << dest[0] << std::endl;
// 验证结果
    float expectedSum = static_cast<float>(N);
    if (dest[0] == expectedSum)
    {
        std::cout << "Result is correct." << std::endl;
    }
    else
    {
        std::cout << "Result is incorrect." << std::endl;
    }
    // 释放内存
    cudaFree(d_src);
    cudaFree(d_dest);
    delete[] src;
    delete[] dest;

    return 0;
}

 1.2.2Work-group内归约优化

Wavefront优化和局部内存优化

由图1-1可知，Naïve Reduction执行时wavefront内部线程存在条件分支，而且对LDS的bank利用率低。首先，Naïve Reduction执行归约的线程ID并不连续，意味着同一个wavefront的线程在kernel执行过程存在条件分支，一部分线程负责归约操作，一部分线程则处于空闲状态。其次，从图1-1的第一层归约可以看出，由于存在空转线程，因而部分bank同样处于空闲状态，LDS的利用率低。

图3-2为改进后的归约算法示意图。如图3-2所示，wavefront内线程不存在条件分支，一个wavefront所能处理的数据将会翻倍，提升wavefront的工作效率，有效地减少了实际工作的wavefronts数目，约为Naïve Reduction的一半。对于局部内存的访问，通过连续的线程访问连续的数据，连续的32个线程将会访问连续的bank，在提升LDS的利用率同时，也有效地避免bank conflict，进一步提升算法性能。完成wavefront优化和局部内存优化的算法版本定义为Divergence-Free Kernel，相对于Naïve Kernel取得良好的性能提升。

图 3-2 Divergence-Free Kernel 归约过程

Fig.3-2 Implementation of Divergence-Free Kernel

循环展开

节针对work-group内归约进行循环展开优化。首先从硬件资源组织上分析，每一个wavefront由64个线程组成(warp由32个线程组成)，wavefront是GPU调度与执行的基本单位，wavefront内所有线程均执行相同的指令，由此可知，在work-group内归约中的for循环中，当运行线程数小于或等于64时，即运行线程都属于同一个wavefront时，可以省去显式的本地同步操作以提升算法性能。

考虑到本文设定work-group内部线程数为256，因此可对for循环进行完全展开，这里需要注意的是，当work-group内实际工作线程的数目大于64(32， NVIDIA GPU)时，仍需要显式的本地同步。

因此，在Divergence-Free Kernel的基础上提出循环展开优化后的work-group内归约优化算法版本Completely-Unroll Kernel，其work-group内归约的伪代码如下：

1.2.3Work-group间归约优化

 归约算法中的work-group间归约主要负责完成对每一个work-group在第二层中得到的局部归约结果的再归约操作，最终得出原始数据集的最终归约结果。work-group间归约总共有三种方法：1)将所有work-group得到的局部归约结果写入到位于global memory中的临时数组中，然后再重新启动归约kernel，进行递归归约操作，直至得到最终归约结果。然而，考虑到启动kernel是一个十分耗时的操作，因此不建议使用。2)将局部归约结果临时数据回传至CPU内存中，在CPU端完成最后的归约操作。但由于数据的回传需要经过PCI-E总线，非常耗时，因此这种方法需要考虑到适当限制开启work-group的数目。3)采用原子操作求得最终的归约结果，这也是本文采用的方法。

本文在归约算法第三层的work-group间归约采用原子操作，主要原因有两点：1) 由于本文实现采用了线程内归约优化，可大大减少开启的work-group数目，从而减少需要进行work-group间归约的局部归约结果数量。2)虽然开启的work-group数目较多，但在GPU目前的调度机制中，能够同时进行work-group间归约，调用原子操作的work-group数量最多为硬件 CU的个数。同时，即使这些work-group在最终执行上，也会存在一定时间间隔，调用原子操作对性能的影响会进一步较小。因此，相对于前两种方法，使用原子操作来完成work-group间归约过程可大大提升整体性能。

CUDA编程-ReduceSum优化记录（文字+图解） - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/628732347

深入浅出GPU优化系列：reduce优化 - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/426978026