數組求和的快速方法（利用cuda的共享內存）--第一部分之源碼分析

代碼來自於這裡

https://code.google.com/p/stanford-cs193g-sp2010/source/browse/trunk/tutorials/sum_reduction.cu

貌似是斯坦福的課程。

核函數分析：

// this kernel computes, per-block, the sum
// of a block-sized portion of the input
// using a block-wide reduction
template
__global__ void block_sum(const DType *input,
		DType *per_block_results,
						const size_t n)
{
	extern __shared__ DType sdata[];

	unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

	// load input into __shared__ memory //一個線程負責把一個元素從全局內存載入到共享內存
	DType x = 0;
	if(i < n)
	{
		x = input[i];
	}
	sdata[threadIdx.x] = x;
	__syncthreads();//等待所有線程把自己負責的元素載入到共享內存

	// contiguous range pattern//塊內進行合並操作，每次合並變為一半.注意threadIdx.x是塊內的偏移，上面算出的i是全局的偏移。
	for(int offset = blockDim.x / 2;
			offset > 0;
			offset >>= 1)
	{
		if(threadIdx.x < offset)//控制只有某些線程才進行操作。
		{
			// add a partial sum upstream to our own
			sdata[threadIdx.x] += sdata[threadIdx.x + offset];
		}

		// wait until all threads in the block have
		// updated their partial sums
		__syncthreads();
	}

	// thread 0 writes the final result//每個塊的線程0負責存放塊內求和的結果
	if(threadIdx.x == 0)
	{
		per_block_results[blockIdx.x] = sdata[0];
	}
}

使用分析：

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

		const size_t block_size = 512;//線程塊的大小。目前有些gpu的線程塊最大為512，有些為1024.
		const size_t num_blocks = (num_elements/block_size) + ((num_elements%block_size) ? 1 : 0);

	// allocate space to hold one partial sum per block, plus one additional
	// slot to store the total sum
		double *d_partial_sums_and_total = 0;//一個線程塊一個和，另外加一個元素，存放所有線程塊的和。
		cudaMalloc((void**)&d_partial_sums_and_total, sizeof(double) * (num_blocks + 1));

		// launch one kernel to compute, per-block, a partial sum//把每個線程塊的和求出來
		block_sum<<>>(d_input, d_partial_sums_and_total, num_elements);

		// launch a single block to compute the sum of the partial sums
	        //再次用一個線程塊把上一步的結果求和。
		//注意這裡有個限制，上一步線程塊的數量，必須不大於一個線程塊線程的最大數量，因為這一步得把上一步的結果放在一個線程塊操作。
		//即num_blocks不能大於線程塊的最大線程數量。
		block_sum<<<1,num_blocks,num_blocks * sizeof(double)>>>(d_partial_sums_and_total, d_partial_sums_and_total + num_blocks, num_blocks);

		// copy the result back to the host
		double device_result = 0;
		cudaMemcpy(&device_result, d_partial_sums_and_total + num_blocks, sizeof(double), cudaMemcpyDeviceToHost);

		std::cout << "Device sum: " << device_result << std::endl;

		// deallocate device memory
		cudaFree(d_input);
		cudaFree(d_partial_sums_and_total);

上面提到這個程序有個限制，即num_blocks不能大於線程塊的最大線程數量。本來這兩者沒有必然聯系的，只是由於程序的設計導致的問題。

這也導致了數組的元素數量也有限制了，不能大於線程塊最大線程數的平方。如果maxThreadsPerBlock=512，那數組數量必須少於512*512；

如果maxThreadsPerBlock=1024，則數組數量必須少於1024*1024。

其實是可以繼續完善這一點的。只需要繼續迭代，保證最後一次是可以在一個線程塊計算即可。

有時間實現一下。