12.2_两遍归约
12.2 两遍归约
该算法包含两个阶段。一个内核执行NumBlocks个并行归约，其中NumBlocks是指用来调用内核的线程块数，得到的结果写入一个中间数组。最终的结果通过一个线程块第二遍调用同一个内核在中间数组基础上生成。代码清单12-1给出了两遍归约内核，它计算一个整数数组的总和。
代码清单12- 1 两遍归约内核
global _void Reduction1_kernel( int *out, const int *in, size_t N) { extern __shared__int sPartials[]; int sum = 0; const int tid = threadIdx.x; for ( size_t i = blockIdx.x*blockDim.x + tid; i < N; i += blockDim.x*gridDim.x ) { sum += in[i]; } sPartials[tid] = sum; _syncthreads(); for ( int activeThreads = blockDim.x>>1; activeThreads; activeThreads >= 1 ) { if ( tid < activeThreads ) { sPartials[tid] += sPartials[tid+activeThreads]; } _syncthreads(); } if ( tid == 0 ) { out[blockIdx.x] = sPartials[0]; } } void Reduction1(int *answer, int *partial, const int *in, size_t N, int numBlocks, int numThreads)unsigned intsharesize $\equiv$ numThreads\*sizeof(int); Reduction1_kernel<< numBlocks，numThreads，shareSize>>>( partial,in,N）; Reduction1_kernel<< 1，numThreads，shareSize $> > >$ { answer，partial，numBlocks）;位于共享内存的数组用来累计每个线程块内的归约。其大小取决于块的线程数，所以当内核启动时必须指定该值。注意：块的线程数必须是2的幂次！
第一个for循环对输入数组落入每个线程中的元素进行求和。如果输入指针被恰当的对齐，由这段代码发起的全部内存事务将被合并，这将最大限度地提高内存带宽。然后，每个线程把它得到的累计值写入共享内存，并在执行对数步长的归约前进行同步操作。
第二个for循环针对共享内存中的值执行对数步长的归约操作。共享内存中前半部分的值被添加到后半部分的值上，而参与的线程数依次减半，直到shared_sum[0]中的一个值包含该块的输出。内核的这部分调用需要该块的线程数是2的幂次。
最后，线程块的输出值写入全局内存。该内核将被调用两次，如主机函数所示：第一次使用N个块，其中N选成可以在输入数组上进行最佳性能的归约操作，然后采用1个块调用内核来累计出最终输出。代码清单12-1展示了调用Reduction1_kernel()的主机函数。需要注意的是，这里需要分配一个用于保存部分和的数组并单独传递。还要注
意，由于内核使用一个未确定大小的共享内存数组，所需的共享内存数量必须由内核指定为<<>>>语法的第三个参数。
CUDA SDK针对这个内核讨论了一些优化方案，重点是减少对数步长的归约过程中的条件代码数量。执行对数步长归约的部分for循环（循环的后部分，当线程数为32或更少时），可以用线程束同步代码实现。由于每个线程块中的线程束是按照锁步方式（lockstep）执行每条指令的，当线程块的活动线程数低于硬件线程束的大小32时，无须再调用__syncthreads()内置函数。由此产生的内核，位于reduction2.cu源代码文件，展示在代码清单12-2中。
重要事项 当编写线程束同步的代码时，必须对来自共享内存的指针使用volatile关键词修饰。否则，编译器引入的优化行为可能会改变内存操作的顺序并且代码将无法正确工作。
代码清单12-2 以循环展开的线程束同步作为结束的归约
__global__void
Reduction2_kernel( int *out, const int *in, size_t N )
{
extern __shared__int sPartials[ ];int sum = 0;
const int tid = threadIdx.x;
for ( size_t i = blockIdx.x*blockDim.x + tid;
i < N;
i += blockDim.x*gridDim.x) {
sum += in[i];
}
sPartials[tid] = sum;
__syncthreads();
for ( int activeThreads = blockDim.x>>1;
activeThreads > 32;
activeThreads >= 1) {
if ( tid < activeThreads) {
sPartials[tid] += sPartials[tid+activeThreads];
}
__syncthreads();
}
if ( threadIdx.x < 32) {
volatile int *wsSum = sPartials;
if ( blockDim.x > 32) wsSum[tid] += wsSum[tid + 32];
wsSum[tid] += wsSum[tid + 16];
wsSum[tid] += wsSum[tid + 8];
wsSum[tid] += wsSum[tid + 4];
wsSum[tid] += wsSum[tid + 2];
wsSum[tid] += wsSum[tid + 1];
if ( tid == 0) {
volatile int *wsSum = sPartials;
out[blockIdx.x] = wsSum[0];
}
}通过把线程数变成一个模板参数，线程束同步的优化方案可以更进一步，从而实现对数步长归约的完全展开。代码清单12-3完整地给出了优化的内核。与Mark Harris的归约演示一脉相承，[1]在编译阶段评估的代码以斜体表示。
代码清单12-3 模板化的完全展开的对数步长的归约
template<unsigned int numThreads> _global__void Reduction3_kernel( int \*out, const int \*in, size_t N)
{ extern __shared__int sPartials[]; const unsigned int tid $=$ threadIdx.x; int sum $= 0$ . for ( size_t i $=$ blockIdx.x\*numThreads $^+$ tid; i $< \mathbb{N}$ i $+ =$ numThreads\*gridDim.x } { sum $+ =$ in[i]; } sPartials[tid] $=$ sum; _syncthreads(); if (numThreads $> = 1024$ ) { if (tid < 512){ sPartials[tid] $+ =$ sPartials[tid + 512]; } _syncthreads(); } if (numThreads $> = 512$ ) { if (tid < 256){ sPartials[tid] $+ =$ sPartials[tid + 256]; } syncthreads(); } if (numThreads $> = 256$ ) { if (tid < 128){ sPartials[tid] $+ =$ sPartials[tid + 128]; } _syncthreads(); } if (numThreads $> = 128$ ) { if (tid < 64){ sPartials[tid] $+ =$ sPartials[tid + 64]; } _syncthreads(); } // warp synchronous at the end if (tid < 32）{ volatile int \*wsSum $=$ sPartials; if (numThreads $> = 64$ ){ wsSum[tid] $+ =$ wsSum[tid + 32); } if (numThreads $> = 32$ ){ wsSum[tid] $+ =$ wsSum[tid + 16); } if (numThreads $> = 16$ ){ wsSum[tid] $+ =$ wsSum[tid + 8); } if (numThreads $> = 8$ ){ wsSum[tid] $+ =$ wsSum[tid + 4); } if (numThreads $> = 4$ ){ wsSum[tid] $+ =$ wsSum[tid + 2); } if (numThreads $> = 2$ ){ wsSum[tid] $+ =$ wsSum[tid + 1); } if (tid $= = 0$ { out [blockIdx.x] $=$ wsSum[0]; }为了实例化代码清单12-3中的函数模板，它必须通过一个独立的主机函数显式地调用。代码清单12-4显示了Reduction3_kernel如何被
另一个函数模板调用，并且主机函数使用了一个switch语句以按照不同线程块大小触发该模板。
代码清单12-4 展开的归约的模板实例化
template<unsigned int numThreads>
void
Reduction3_template( int \*answer, int \*partial, const int \*in, size_t N, int numBlocks)
{ Reduction3_kernel<numThreads><\n numBlocks, numThreads, numThreads\*sizeof(int) $\ggg$ ( partial, in, N); Reduction3_kernel<numThreads><\n 1, numThreads, numThreads\*sizeof(int) $\ggg$ ( answer, partial, numBlocks);
}
void
Reduction3( int \*out, int \*partial, const int \*in, size_t N, int numBlocks, int numThreads)
{ switch ( numThreads ) { case 1: return Reduction3_template< 1>(...); case 2: return Reduction3_template< 2>(...); case 4: return Reduction3_template< 4>(...); case 8: return Reduction3_template< 8>(...); case 16: return Reduction3_template< 16>(...); case 32: return Reduction3_template< 32>(...); case 64: return Reduction3_template< 64>(...); case 128: return Reduction3_template< 128>(...); case 256: return Reduction3_template< 256>(...); case 512: return Reduction3_template< 512>(...); case 1024: return Reduction3_template<1024>(...); }