CUDA规约求和/最大值

1. 规约思想

假设有一个复杂的问题P，而它看起来与一个已知的问题Q很相似，可以试着在两个问题间找到一个归约(reduction, 或者transformation)。

对于问题的先后，归约可以达到两个目标：
(1)已知Q的算法,那么就可以把使用了Q的黑盒的P的解决方法转化成一个P的算法。
(2)如果P是一个已知的难题，或者特别地,如果P的下限，那么同样的下限也可能适用于Q。前一个归约是用于获取P的信息;而后者则是用于获取Q的信息
归约让我们理解两个问题间的关系，它是一种技术，用于寻找解决某个问题或它的变形。

2. 规约求和/最大值

规约求最大值和规约求和的思想一样，以下将通过规约求和的思想进行阐述。

规约求和（Reduction Summation）是一种计算方法，通常用于向量或数组求和。它通过将数据划分为更小的块，然后对每个块进行求和，最后将所有块的求和结果相加得到最终结果。

具体实现过程可能因语言和环境而异，但通常包括以下步骤：

1. 将输入数据划分为更小的块。
2. 对每个数据块进行求和。
3. 将所有块的求和结果相加得到最终结果。

在规约求和中，加法的交换律和结合律可以保证求和的正确性，无论数据的排列顺序如何，只要保证相同元素在同一位置上，就可以得到正确的结果。

下面我们举一个简单的例子来说明规约思想，假设数组X有8个元素，现要使用规约思想求其元素和。每个步骤都将数据两两分组，然后并行计算每组的元素和，最后得到一个结果，这就是规约的过程：

3. CUDA实现规约求和/求最大值

CUDA是为并行计算而生的，使用CUDA可以很容实现上述的数组规约求和算法。不过有一点需要注意，就是必须确保每个步骤的所有线程是同步的，也即所有线程计算完成之后再进入下一步骤的计算，否则会导致结果错误。
在CUDA中，可以调用__syncthreads函数方便地同步同一个线程块中的所有线程，因此我们可以使用同一个线程块中的多个线程做规约运算。那么问题来了，如果数据量很大，一个线程块不能完成所有数据的规约运算该怎么办呢？答案是分块处理，将数据平均分成多个部分，每部分都分配给一个线程块做规约运算。因此每个线程块最后得到一个规约结果，最后再将多个规约结果求和/求最大值，即得到最后结果。如下图所示：

值得注意的是，最后的多个规约结果可以拷贝到CPU中计算，也可以在GPU进行计算。在GPU中进行计算时，当规约结果数量小于单个线程块的线程数时，即可得到最终结果。

4. 编程实现

编程思路：首先把输入数组划分为更小的数据块，之后用一个线程计算一个数据块的部分和，最后把所有部分和再求和得出最终结果。

上面演示了单个线程块内的归约求和过程，结合下面的代码解释：

第一轮循环，对应上图的第二行，数组中索引为偶数的值被对应部分和替代。
第二轮循环，对应上图第三行，数组中索引为 4 的整数倍的值被对应部分和替代。
第三轮循环，对应上图第四行，数组中索引为 8 的整数倍的值被对应部分和替代。
…
直至 stride 的值为 blockDim.x 的一半，这一步块内和就计算完了，并存储在threadIdx.x 为 0 的线程中。但对于不同块内的和，就需要将每个块内 threadIdx.x 为 0 的线程的和相加才能求得总的向量和。

规约求最大值核函数

__global__ void reduceMax(float* maxGpu, float *d_A, int n)
{
	unsigned int t = threadIdx.x;
	__shared__ float partialSum[THREAD_LENGTH];
	if (blockIdx.x*blockDim.x + t < n)
		partialSum[t] = d_A[blockIdx.x*blockDim.x + t];
	else
		partialSum[t] = 0;
	__syncthreads();  //将数组加载到共享存储器。
	for (unsigned int stride = 1; stride < blockDim.x; stride *= 2)
	{
		if (t % (2 * stride) == 0)
		{
			//partialSum[t] += partialSum[t + stride]; //规约求和
			if (partialSum[t] < partialSum[t + stride])//规约求最大值
				partialSum[t] = partialSum[t + stride];
		}
		__syncthreads();
	}
	/*将每个线程块内 threadIdx.x 为零的线程中的值传回d_A。
	主机函数中将会对这几个线程求和，以得到最终的和。*/
	if (t == 0) {
		//printf("blockIdx x=%d\n", blockIdx.x);
		maxGpu[blockIdx.x] = partialSum[t];
	}
}

全程在GPU中计算最大值

    ImagePro* tempGpu = nullptr;
	cudaMalloc((void **)&tempGpu, w*h * sizeof(ImagePro));
	cudaMemcpy(tempGpu, pSrc, w*h * sizeof(ImagePro), cudaMemcpyDeviceToDevice);

	int n = w*h;
	int blocks = ceil((double)n / THREAD_LENGTH);
	ImagePro *maxArrayGpu = nullptr;
	ImagePro *minArrayGpu = nullptr;
	cudaMalloc((void **)&maxArrayGpu, blocks * sizeof(ImagePro));
	cudaMalloc((void **)&minArrayGpu, blocks * sizeof(ImagePro));
	while (n > THREAD_LENGTH) {
		blocks = ceil((double)n / THREAD_LENGTH);
		//printf("blocks=%d\n", blocks);
		reduceMax << <blocks, THREAD_LENGTH >> > (maxArrayGpu, tempGpu, n);
		cudaDeviceSynchronize();
		cudaMemcpy(tempGpu, maxArrayGpu, blocks * sizeof(ImagePro), cudaMemcpyDeviceToDevice);
		cudaDeviceSynchronize();
		n = blocks;
	}
	blocks = ceil((double)n / THREAD_LENGTH);
	//printf("blocks=%d\n", blocks);
	reduceMax << <blocks, THREAD_LENGTH >> > (maxArrayGpu, tempGpu, n);
	cudaDeviceSynchronize();
	//把结果拷贝到CPU进行check
	cudaMemcpy(maxValueCpu, maxArrayGpu, 1 * sizeof(ImagePro), cudaMemcpyDeviceToHost);
	cudaDeviceSynchronize();

CPU完成最后的计算

    reduceMax << <ceil((double)n / THREAD_LENGTH), THREAD_LENGTH >> >(maxArrayGpu, tempGpu, n);
	LAST_KERNEL_CHECK();
	cudaMemcpy(maxValueCpu, maxArrayGpu, n * sizeof(ImagePro), cudaMemcpyDeviceToHost); 
	LAST_KERNEL_CHECK();
	for (int i = 0; i < ceil((double)n / THREAD_LENGTH); ++i)
	{  
		if (maxValue < maxValueCpu[i*THREAD_LENGTH])
			maxValue = maxValueCpu[i*THREAD_LENGTH];
	}

5. 参考文献

[1] CUDA----规约
[2] CUDA 归约求和算法（一）
[3] CUDA加速——基于规约思想的数组元素求和
[4] CUDA编程 | 规约算法求最值、求和