让我们分块看代码,一路回答你的问题:
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];
}
上述代码遍历大小为N 的数据集。为了便于理解,我们可以做出的假设是N > blockDim.x*gridDim.x,最后一项就是网格中的线程总数。由于N 大于总线程数,因此每个线程都在对数据集中的多个元素求和。从给定线程的角度来看,它是对由线程的网格维度分隔的元素求和 (blockDim.x*gridDim.x) 每个线程将其总和存储在名为 sum 的本地(可能是寄存器)变量中。
sPartials[tid] = sum;
__syncthreads();
随着每个线程完成(即,因为它的 for 循环超过 N),它会将其中间的 sum 存储在共享内存中,然后等待块中的所有其他线程完成。
for ( int activeThreads = blockDim.x>>1;
activeThreads > 32;
activeThreads >>= 1 ) {
if ( tid < activeThreads ) {
sPartials[tid] += sPartials[tid+activeThreads];
}
__syncthreads();
}
到目前为止,我们还没有讨论过方块的尺寸——这无关紧要。假设每个块都有 32 个线程的整数倍。下一步将开始将存储在共享内存中的各种中间和收集到越来越小的变量组中。上面的代码首先选择线程块(blockDim.x>>1)中的一半线程,并使用这些线程中的每一个来组合共享内存中的两个部分和。因此,如果我们的线程块从 128 个线程开始,我们只需使用其中的 64 个线程将 128 个部分和减少为 64 个部分和。这个过程在 for 循环中重复地继续,每次将线程切成两半并合并部分和,每个线程一次两个。只要activeThreads > 32,这个过程就会继续。所以如果activeThreads 是64,那么这64 个线程会将128 个部分和组合成64 个部分和。但是当activeThreads 变为 32 时,for 循环终止,没有将 64 个部分和合并为 32。所以在完成这段代码时,我们取了(32 个线程的任意倍数) ) 线程块,并将我们开始时的许多部分和减少到 64。这个将 256 个部分和、128 个部分和、64 个部分和组合在一起的过程必须在每次迭代时等待所有线程(在多个 warp 中)完成他们的工作,所以__syncthreads(); 语句在for循环的每遍中执行。
请记住,此时,我们已将线程块减少到 64 个部分和。
if ( threadIdx.x < 32 ) {
对于此后的内核剩余部分,我们将只使用前 32 个线程(即第一个 warp)。所有其他线程将保持空闲。请注意,在此之后也没有__syncthreads();,因为这将违反使用它的规则(所有线程都必须参与__syncthreads();)。
volatile int *wsSum = sPartials;
我们现在正在创建一个指向共享内存的volatile 指针。理论上,这告诉编译器它不应该进行各种优化,例如将特定值优化到寄存器中。为什么我们以前不需要这个?因为__syncthreads(); 也带有a memory-fencing function。 __syncthreads(); 调用除了导致所有线程在屏障处相互等待之外,还会强制所有线程更新回共享内存或全局内存。然而,我们不能再依赖这个特性了,因为从现在开始我们将不再使用__syncthreads();,因为我们已经将自己(对于内核的其余部分)限制为单个 warp。
if ( blockDim.x > 32 ) wsSum[tid] += wsSum[tid + 32]; // why do we need this
之前的归约块给我们留下了 64 个部分和。但我们此时将自己限制为 32 个线程。因此,我们必须再进行一次组合,将 64 个部分和合并为 32 个部分和,然后才能继续进行剩余的归约。
wsSum[tid] += wsSum[tid + 16]; //how these statements are executed in paralle within a warp
现在我们终于进入了一些 warp 同步编程。这行代码取决于 32 个线程同步执行的事实。为了理解为什么(以及它是如何工作的),将其分解为完成这行代码所需的操作序列会很方便。它看起来像:
read the partial sum of my thread into a register
read the partial sum of the thread that is 16 higher than my thread, into a register
add the two partial sums
store the result back into the partial sum corresponding to my thread
所有 32 个线程都将按照上述顺序锁步。所有 32 个线程都将从将wsSum[tid] 读入(线程本地)寄存器开始。这意味着线程 0 读取 wsSum[0],线程 1 读取 wsSum[1] 等。之后,每个线程将 another 部分和读入不同的寄存器:线程 0 读取 wsSum[16],线程 1 读取 @987654355 @ 等。确实,我们并不关心wsSum[32](及更高)的值;我们已经将它们折叠成前 32 个 wsSum[] 值。然而,正如我们将看到的,只有前 16 个线程(在这一步)对最终结果有贡献,所以前 16 个线程将把 32 个部分和合并为 16。接下来的 16 个线程也将起作用,但他们只是在做垃圾工作——它将被忽略。
上述步骤将 32 个部分和合并到 wsSum[] 中的前 16 个位置。下一行代码:
wsSum[tid] += wsSum[tid + 8];
以 8 的粒度重复这个过程。同样,所有 32 个线程都处于活动状态,微序列是这样的:
read the partial sum of my thread into a register
read the partial sum of the thread that is 8 higher than my thread, into a register
add the two partial sums
store the result back into the partial sum corresponding to my thread
所以前 8 个线程将前 16 个部分和 (wsSum[0..15]) 组合成 8 个部分和(包含在 wsSum[0..7] 中)。接下来的 8 个线程也将 wsSum[8..23] 合并到 wsSums[8..15],但是写入 8..15 发生在线程 0..8 读取这些值之后,因此有效数据不会损坏。这只是额外的垃圾工作。对于经线中的其他 8 个线程块也是如此。因此,此时我们已将部分兴趣总和合并到 8 个位置。
wsSum[tid] += wsSum[tid + 4]; //this combines partial sums of interest into 4 locations
wsSum[tid] += wsSum[tid + 2]; //this combines partial sums of interest into 2 locations
wsSum[tid] += wsSum[tid + 1]; //this combines partial sums of interest into 1 location
而这几行代码遵循与前两行类似的模式,将 warp 分成 8 组,每组 4 个线程(只有第一个 4 线程组对最终结果有贡献),然后将 warp 分成 16 组,每组 2 个线程,只有第一个 2 线程组对最终结果有贡献。最后,分成 32 组,每组 1 个线程,每个线程生成一个部分和,只有第一个部分和是感兴趣的。
if ( tid == 0 ) {
volatile int *wsSum = sPartials;// why this statement is needed?
out[blockIdx.x] = wsSum[0];
}
最后,在上一步中,我们将所有部分和缩减为一个值。现在是时候将单个值写入全局内存了。我们完成了减少吗?也许,但可能不是。如果上述内核仅使用 1 个线程块启动,那么我们就完成了——我们最终的“部分”总和实际上是数据集中所有元素的总和。但是如果我们启动了多个区块,那么每个区块的最终结果仍然是一个“部分”和,并且所有区块的结果必须相加(以某种方式)。
然后回答你的最后一个问题?
我不知道为什么需要那个声明。
我的猜测是它是在缩减内核的先前迭代中遗留下来的,程序员忘记删除它,或者没有注意到它不需要它。也许其他人会知道这个问题的答案。
最后,cuda reduction sample 提供了非常好的参考代码供学习,随附的pdf document 很好地描述了可以在此过程中进行的优化。