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();