代码来自于这里
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<class DType>
__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<<<num_blocks,block_size,block_size * sizeof(double)>>>(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。
其实是可以继续完善这一点的。只需要继续迭代,保证最后一次是可以在一个线程块计算即可。
有时间实现一下。