4.7 优化应用举例
4.7.2 并行归约的优化
本节将通过 SDK 中的 reduction 示例演示对 CUDA 程序进行逐步优化的过程。reduction 并行归约（或叫并行缩减）提供了一种可以实现并行求和、求乘、求最大（小）值等各种操作 的方法，用途十分广泛。本节的 reduction 示例中使用的运算是加，可以实现对数组元素的并 行求和操作。下面给出 7 个版本的 kernel 程序，并逐一进行分析。在这一段代码中，假设有 6 个一维 block，每个 block 有 512 个线程，而输入数组中的元素数量为 3072 个，由每个线程对 应处理一个数据。
4.7.2.1 reduce0
以下是 reduce0 版本的 kernel 程序。
//---reduction_kernel.cu---
#ifndef _REDUCE_KERNEL_H_
#define _REDUCE_KERNEL_H_
#include <stdio.h>
#include "sharedmem.cuh"
#ifdef __DEVICE_EMULATION__
#define EMUSYNC __syncthreads()
#else
#define EMUSYNC
#endif
// Macros to append an SM version identifier to a function name
// This allows us to compile a file multiple times for different architecture // versions
// The second macro is necessary to evaluate the value of the SMVERSION macro // rather than appending "SMVERSION" itself
#define FUNCVERSION(x, y) x ## _ ## y
#define XFUNCVERSION(x, y) FUNCVERSION(x, y)
#define FUNC(NAME) XFUNCVERSION(NAME, SMVERSION)
/*
template <class T>
__global__ void
FUNC(reduce0)(T *g_idata, T *g_odata) {
__syncthreads();
}
//将 block 的计算结果写回全局存储器 if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
上段代码中,首先由每个线程读入一个数据，完成从全局存储器向共享存储器的数据拷贝，
随后进行了一次同步操作，以保证所有数据都可以安全地被其他线程访问。在接下来的 for() 循环中，每一轮都只使用上一轮循环中一半的线程进行 tid 与 tid+s 的求和，即规约的过程。
在这里，s 可以理解为跨度，对于每一个 block，有：
第一次循环（s=1），只有 tid=0、2、4、6、8…510 线程真正执行计算，即 sdata[tid]
与其后跨度为 1 的元素的加和操作。
第二次循环（s=2），只有 tid=0、4、8、12、16…508 线程真正执行计算，即 sdata[tid]
与其后跨度为 2 的元素的加和。
依此类推。
直到最后一次（s=256），只有 tid=0 线程真正执行计算，即 sdata[tid]与其后跨度为 256 的元素的加和。此时 sdata[0]中的结果即为该 block 输入的 512 个数据之和。
注意到，每次循环中都进行了一次__syncthreads()同步，这是因为下一轮循环中的线程要 用到上一轮中其他线程计算的结果。reduce0 的执行过程如图 4-17 所示。
图 4-17 reduce0 加和过程 这一版本代码中，主要存在两个问题：
（1）在判断线程是否需要参与运算时使用了取模运算，由于 GPU 的整数处理单元功能 较弱，做整数的模运算和除法运算的开销都特别大，应尽量避免使用或用位运算代替。
（2）每次循环时，只有 tid 是 2*s 倍数的线程才真正参与运算，因此每个 warp 中都只有 部分线程执行加和操作（s=1 时 warp 中参与运算的线程最多，但也只有一半），但同时绝大多 数 warp 都有线程要执行加和。于是，在分支时，绝大多数的 warp 都有要执行加法的线程，需 要执行加法分支，只有少数的 warp 中的 thread 全部为空分支。
4.7.2.2 reduce1
针对 reduce0 中的“取模”问题，我们给出 reduce1 版本。
/* 这个版本使用连续线程，但它交叠的取址方式导致严重的分区冲突问题 */
template <class T>
__global__ void
FUNC(reduce1)(T *g_idata, T *g_odata) {
同 reduce0
for(unsigned int s=1; s < blockDim.x; s *= 2) {
int index = 2 * s * tid;
if (index < blockDim.x) {
sdata[index] += sdata[index + s];
}
__syncthreads();
}
//将本 block 的计算结果写回全局存储器 if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}
下面对 for()循环进行分析。对于每一个 block：
第一次循环（s=1），index<blockDim.x，即要求 2*tid<512，那么只有 tid=0、2、4、6…254 线程在做有效计算，sdata[2*tid]与其后跨度为 1 的元素加和。
第二次循环（s=2），4*tid<512，那么只有 tid=0、4、8、12…128 在做有效计算，sdata[4*tid]
与其后跨度为 2 的元素加和。
依此类推。
直至最后一次循环（s=256），那么只有 tid=0 在做有效运算，sdata[0]+=sdata[256]。
其实，reduce0 与 reduce1 的计算逻辑是一样的，reduce1 通过引入 index 避免了取模运算。
但是在 reduce0、reduce1 两个版本中，随着 s 的不断增加，都会造成越来越剧烈的 bank conflict。
例如 T 为 float 类型，那么 s=1 时，0～15 号线程 half-warp 分别访问 sdata[1/3/5/../31]，形成 2 路冲突，图 4-18 中下划线标注；s=2 时，half-warp 会形成 4 路冲突，图中以带阴影的矩形进 行标注；s=4 时，half-warp 会形成 8 路冲突，图中以不带阴影的矩形标注；以此类推。
图 4-18 reduce1 共享存储器 bank conflict 示意
4.7.2.3 reduce2
针对前两个版本中存在的 bank conflict 及 warp 分支问题，我们进一步分析 reduce2 版本。
/*
这个版本的 reduction 使用序列地址，避免了 bank conflict.
*/
template <class T>
__global__ void
FUNC(reduce2)(T *g_idata, T *g_odata) {
reduce2 版本中，每个 thread 操作的元素都是相邻的，因此不会造成 bank conflict。例如对 于每一个 block，第一次循环（s=256），tid=0…255 都会进行 sdata[tid]与 sdata[tid+256]元素的 加和，所以说相邻 thread 操作相邻元素。此外，每个 warp 中的所有线程要么全部执行加法，
要么全部不执行，也避免了分支效率问题。
4.7.2.4 reduce3
下面的 reduce3 版本使得相同的线程数可以完成更多元素的求和操作。上面的第三个版本 已经将 shared memory 内的规约算法的效率挖掘殆尽了。剩下的工作就是在编程技巧和 shared memory 外做文章了。注意到上面的算法中，每次循环中需要访问的 shared memory 和参与运 算的线程都比上一次要少。
下面的 reduce3 版本在访问 global memory 时就进行了一次加法，这样每个 block 只需要同 样的线程，就能完成多一倍的元素的规约求和操作。
/*
这个版本的程序使用了 n/2 个线程--This version uses n/2 threads -- 在从全局内存读的时候执行了第一层次的归约
*/
template <class T>
__global__ void
FUNC(reduce3)(T *g_idata, T *g_odata) {
SharedMemory<T> smem;
T *sdata = smem.getPointer();
//执行第一层次的归约（并行读数据），写入共享存储器 个数据（而不是原来的 512 个），例如对于线程 0 有 sdata[0]=g_idata[0]+g_idata[511]；第 1 个 block 负责 1024～2047 的数据加和；第 5 个 block 负责 5120～6144 的数据加和，而 6144 是之 前能处理的数据即 3072 的两倍。
4.7.2.5 reduce4
根据 warp 内线程不需要同步的特点，我们进一步得到 reduce4 版本。这个版本中每次循 环都进行了一次__syncthreads()同步。实际上，当只剩最后 32 个线程时，也就只有一个 warp 需要执行加法，其他的 warp 都执行空分支。由于一个 warp 内的运算总是满足顺序一致性的，
因此在一个 warp 中也就不需要再进行同步了。
注意到条件编译#ifdef __DEVICE_EMULATION，这是为了减少 CPU 模拟运行时的分支。
此外，由于 emu 模式下 CPU 的 warp size 是 1，为了能够得到正确的模拟结果，仍然需要在最 后一个 warp 中加入 EMUSYNC。
/*
这个版本的代码展开了最后一个 warp 以避免同步操作
*/
template <class T, unsigned int blockSize>
__global__ void
FUNC(reduce4)(T *g_idata, T *g_odata) {
同 reduce3
for(unsigned int s=blockDim.x/2; s>32; s>>=1) {
if (tid < s)
{
循环在 GPU 上的执行效率不高，reduce5 版本完全展开了循环。当每个 block 中的 thread 数量是 2 的 n 次方时，下面的代码可以获得最高的效率。
/*
这个版本的代码完全展开了循环。使用了一个模板参数进行代码优化（线程数需是 2 的幂次方）。这 要求主机端有一个 switch 语句来处理不同大小的线程块。
*/
template <class T, unsigned int blockSize>
__global__ void
FUNC(reduce5)(T *g_idata, T *g_odata) {
最后的版本中在循环读取 global memory，并将读到的数据直接加到 shared memory 中。
在这个循环中，绝大多数 warp 都能保持满载，因此可以提高效率。
/*
这个版本的代码中，每个线程顺序地加若干个元素，这就减少了整个算法的开销，而整个复杂度还维持在 O(n)（step complexity 为 O(logn)）
*/
template <class T, unsigned int blockSize>
__global__ void
FUNC(reduce6)(T *g_idata, T *g_odata, unsigned int n) {
template <class T>
void
FUNC(reduce)(int size, int threads, int blocks,
int whichKernel, T *d_idata, T *d_odata)
case 4:
switch (threads) {
case 512:
FUNC(reduce4)<T, 512><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
略//256、128、…、2 case 1:
FUNC(reduce4)<T, 1><<< dimGrid, dimBlock, smemSize >>>(d_idata, d_odata); break;
} break;
略//case 5、6、default 与 case4 类同 }
extern "C"
void FUNC(reduceInt)(int size, int threads, int blocks,
int whichKernel, int *d_idata, int *d_odata) {
FUNC(reduce)<int>(size, threads, blocks, whichKernel, d_idata, d_odata);
}
//同理，有 FUNC(reduceFloat)、FUNC(reduceDouble)等
#endif