数组求和的快速方法（利用cuda的共享内存）--第二部分之程序完善

上一篇提到，那份源码的使用是有限制的。

这次来完善一下。其实就是迭代多次，使得最后一次刚好在一个线程块可以求和。

完善部分：

template<class DType>
DType array_sum_gpu(DType *dev_array,const int array_size,DType *dev_result)
{
	//const size_t max_block_size  = 512;//目前有些gpu的线程块最大为512，有些为1024.
	const size_t block_size = 512;//线程块的大小。

	size_t num_elements = array_size;
	size_t num_blocks = (num_elements/block_size) + ((num_elements%block_size) ? 1 : 0);

	double *dev_input_array = 0;
	double *dev_block_sums = 0;//一个线程块一个和。


	while(num_elements > block_size)
	{
		if(dev_block_sums == 0)//第一次
		{
			dev_input_array = dev_array;
		}
		else //除了第一次
		{
			if(dev_input_array != dev_array)
				cudaFree(dev_input_array);
			dev_input_array = dev_block_sums;
		}
		num_blocks = (num_elements/block_size) + ((num_elements%block_size) ? 1 : 0);
		//给输出结果分配内存
		cudaMalloc((void**)&dev_block_sums, sizeof(double) * (num_blocks ));
		// launch one kernel to compute, per-block, a partial sum//把每个线程块的和求出来
		block_sum<<<num_blocks,block_size,block_size * sizeof(double)>>>(dev_input_array, dev_block_sums, num_elements);
		num_elements = num_blocks;
	}

	block_sum<<<1,num_elements,num_elements * sizeof(double)>>>(dev_block_sums, dev_result, num_elements);
	double result = 0;
	cudaMemcpy(&result, dev_result, sizeof(double), cudaMemcpyDeviceToHost);
	cudaFree(dev_block_sums);
	return result;

}

核函数block_sum还是原来的代码。

下面是测试我的代码;

void test_sum2()
{
	// create array of 256k elements
	//const int num_elements = 1<<18;//=512*512=262144
	const int num_elements = 1<<20;

	// generate random input on the host
	std::vector<double> h_input(num_elements);
	for(int i = 0; i < h_input.size(); ++i)
	{
		h_input[i] = 1;//random_num<double>();
	}

	const double host_result = std::accumulate(h_input.begin(), h_input.end(), 0.0f);
	std::cerr << "Host sum: " << host_result << std::endl;

	// move input to device memory//分配内存
	double *d_input = 0;
	cudaMalloc((void**)&d_input, sizeof(double) * num_elements);
	cudaMemcpy(d_input, &h_input[0], sizeof(double) * num_elements, cudaMemcpyHostToDevice);

	double *dev_result=0;
	cudaMalloc((void**)&dev_result,sizeof(double));

	double sum = array_sum_gpu(d_input,num_elements,dev_result);
	std::cout << "Device sum: " << sum << std::endl;


}

可以把数组数量num_elements调到很大，代码仍然能运行正确。但是如果是之前那份代码，就不可以了。

其实这个程序还是有点限制的。

请注意第一次求num_blocks.

size_t num_blocks = (num_elements/block_size) + ((num_elements%block_size) ? 1 : 0);

万一第一次求出的num_blocks大于线程块的最大数量，一般是65535，那就不行了。

所以如果数组的元素数量大于1024*65535，那就无法计算了。

解决这中问题的通常方法，是让一个线程串行执行多个相同任务。

由于求解的问题暂时没有超过这个数量级（6-7千万），所以先这样。