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: %d\n", 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();