cuda第一次计算耗时_CUDA编程入门（四）并行归约算法

这一篇我们一起学习一下如何使用CUDA实现并行归约算法。

首先我们要知道什么是并行归约。并行归约（Reduction）是一种很基础的并行算法，简单来说，我们有N个输入数据，使用一个符合结合律的二元操作符作用其上，最终生成1个结果。这个二元操作符可以是求和、取最大、取最小、平方、逻辑与或等等。

我们以求和为例，假设输入如下：

int array[8] =  [3, 1, 7, 0, 4, 1, 6, 3]

在串行的情况下，算法很容易实现，一般我们会使用下面这样的代码。

int sum = 0;
for (int i = 0; i < N; i++)
    sum += array[i]

但当我们试图进行并行计算时，问题就会变得复杂。

由于加法的交换律和结合律，数组可以以任意顺序求和。所以我们会自然而然产生这样的思路：首先把输入数组划分为更小的数据块，之后用一个线程计算一个数据块的部分和，最后把所有部分和再求和得出最终结果。

依照以上的思路，我们可以按下图这样计算。就上面的输入例子而言，首先需要我们开辟一个8个int的存储空间，图中一行的数据代表我们开辟的存储空间。计算时首先将相邻的两数相加（也称相邻配对），结果写入第一个数的存储空间内。第二轮迭代时我们再将第一次的结果两两相加得出下一级结果，一直重复这个过程最后直到我们得到最终的结果，空白方框里面存储的内容是我们不需要的。这个过程相当于，每一轮迭代后，选取被加数的跨度翻倍，后面我们核函数就是这样实现的。

相邻配对实现并行归约

相比与串行计算，我们只用了3轮迭代就得出了8个数的和，时间复杂度由O(N)变为O(logN)。

当然使用归约算法时我们不会只有这么小的输入数组，这时候就要用到线程块了，我们可以把输入数组先划分成很多包含8个int型值的数组，每一个小数组分配给一个线程块，最后再将所有的结果传回主机串行求和。

主函数如下，与我们上一个例程结构类似，先初始化，分配内存，然后运行核函数，最后和CPU对照组对比，检验结果是否正确。我们重点关注分配线程网格和线程块的主干部分。与上面例子里讲的8元素数组不同，实际使用时为了达到高加速比，我们会把输入数组分成更大的块。

这里我们使用一个16M（16777216）大小的输入数组，分成的小数组大小为1K（1024）。所以程序里将block设置为（1024, 1）的大小，每个线程块完成一个小数组的求和。一共需要16K（16384）个线程块，所以我们设置grid为（16384, 1）。每个线程块的求和结果都保存在全局内存里，等所有线程完成后，统一传回主机，在主机里串行求和得到最后的结果。注意这里的block和grid设置并不是最优，只是为了简单，下一篇中我们会进行优化。

int main(int argc,char** argv)
{
    initDevice(0);

    //initialization

    int size = 1 << 24;
    printf("    with array size %d  ", size);

    //execution configuration
    int blocksize = 1024;
    if (argc > 1)
    {
        blocksize = atoi(argv[1]);   //从命令行输入设置block大小
    }
    dim3 block(blocksize, 1);
    dim3 grid((size - 1) / block.x + 1, 1);
    printf("grid %d block %d n", grid.x, block.x);

    //allocate host memory
    size_t bytes = size * sizeof(int);
    int *idata_host = (int*)malloc(bytes);
    int *odata_host = (int*)malloc(grid.x * sizeof(int));
    int * tmp = (int*)malloc(bytes);

    //initialize the array
    initialData_int(idata_host, size);

    memcpy(tmp, idata_host, bytes);
    double timeStart, timeElaps;
    int gpu_sum = 0;

    // device memory
    int * idata_dev = NULL;
    int * odata_dev = NULL;
    CHECK(cudaMalloc((void**)&idata_dev, bytes));
    CHECK(cudaMalloc((void**)&odata_dev, grid.x * sizeof(int)));

    //cpu reduction 对照组
    int cpu_sum = 0;
    timeStart = cpuSecond();
    //cpu_sum = recursiveReduce(tmp, size);
    for (int i = 0; i < size; i++)
        cpu_sum += tmp[i];
    timeElaps = 1000*(cpuSecond() - timeStart);

    printf("cpu sum:%d n", cpu_sum);
    printf("cpu reduction elapsed %lf ms cpu_sum: %dn", timeElaps, cpu_sum);


    //kernel reduceNeighbored

    CHECK(cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice));
    CHECK(cudaDeviceSynchronize());
    timeStart = cpuSecond();
    reduceNeighbored <<<grid, block >>>(idata_dev, odata_dev, size);
    cudaDeviceSynchronize();
    cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
    gpu_sum = 0;
    for (int i = 0; i < grid.x; i++)
        gpu_sum += odata_host[i];   
    timeElaps = 1000*(cpuSecond() - timeStart);

    printf("gpu sum:%d n", gpu_sum);
    printf("gpu reduceNeighbored elapsed %lf ms     <<<grid %d block %d>>>n",
        timeElaps, grid.x, block.x);

    // free host memory

    free(idata_host);
    free(odata_host);
    CHECK(cudaFree(idata_dev));
    CHECK(cudaFree(odata_dev));

    //reset device
    cudaDeviceReset();

    //check the results
    if (gpu_sum == cpu_sum)
    {
        printf("Test success!n");
    }
    return EXIT_SUCCESS;
}

核函数具体编程实现如下。我们利用一个stride变量，实现不同轮数时被加数的选择。每当计算一轮后，选取被加数的跨度会乘2。在这里，有心的同学就会发现了，第一轮计算中，其实只有一半的线程是活跃的，而且每进行一轮计算后，活跃的线程数都会减少一半，这是条件表达式的使用造成的。比如在第一轮迭代时，只有偶数ID的线程会为True，其主体才能得到执行。这会导致线程束的分化，也就是说只有一部分线程是活跃的，但是所有的线程仍然都会被调度，因为硬件调度线程是以线程束（连续32个线程）为单位进行调度的。这肯定会影响我们程序的执行效率，不过不用太担心，我们下一篇就会提出解决这个问题的方法。

__global__ void reduceNeighbored(int * g_idata,int * g_odata,unsigned int n) 
{
    //set thread ID
    unsigned int tid = threadIdx.x;
    //boundary check
    if (tid >= n) return;
    //convert global data pointer to the 
    int *idata = g_idata + blockIdx.x*blockDim.x;
    //in-place reduction in global memory
    for (int stride = 1; stride < blockDim.x; stride *= 2)
    {
        if ((tid % (2 * stride)) == 0)
        {
            idata[tid] += idata[tid + stride];
        }
        //synchronize within block
        __syncthreads();
    }
    //write result for this block to global mem
    if (tid == 0)
        g_odata[blockIdx.x] = idata[0];

}

最后我们看看运行结果，cpu花费89.3ms，而gpu花费2.5ms，运算结果一致。

终端截图

完整代码和编译生成的可执行文件放在这里：

https://github.com/ZihaoZhao/CUDA_study/tree/master/Reduction

总结一下，这一篇我们使用CUDA完成了一个基础的并行归约算法，不过这算法还有很大的优化空间，接下来我们一起对这个算法进行优化，看看能再加快多少。