束内规约与块内规约问题
写在前面：规约问题在 cuda 编程中应用非常广泛，笔者最近在研究 faster transformer 源码，趁此机会结合 nivida 官方的代码对规约手段进行总结。
1 应用背景 关于规约的定义，相信能读到这篇文章的读者都不陌生，笔者在早期的文章中也介绍过一些规约方法，基本思想都是折半规约，主要应用于较大元素规模的向量规约，有兴趣的读者可以移步【cuda编程】cuda编程中的并行规约问题。
本文要介绍的规约场景与之前有所不同，主要应用于矩阵规约，也就是说本文假设的输入变量的维度是 2 维的，形状为 [batch_size, hidden_units]，规约之后的输出变量形状为 [batch_size, ]。
接下来，本文将以规约求和为例介绍两种规约方式：束内规约、块内规约。
2 束内规约 束内规约，也就是在一个线程束内对某个变量进行规约。我们知道 cuda 架构下指令是以线程束（相邻的 32 个线程）为基本单元执行的，线程束内也可以通过束内洗牌指令进行通信，所以这提供了一个很好的束内规约思路。下面是 nvidia 提供的基础的一个规约设备函数。
template __inline__ __device__t warpreducesum(t val){ for(int mask = 16; mask > 0; mask >>= 1) val += __shfl_xor_sync(final_mask, val, mask, 32); return val;} 这个设备函数可以求出当前线程所在线程束的指定变量的规约和，原理涉及洗牌指令的计算逻辑，不再赘述。
当矩阵宽度 hidden_units 较小时，通常可以使用一个 warp 处理一行数据，一个 block 内可以处理多行数据，笔者给出具体的核函数如下：
// 一个 warp 处理一行数据template__global__ void matrix2dwarpreducesum(const t* inp, t*out, const uint32_t hidden_units) { uint32_t tid = threadidx.x; uint32_t lane_id = tid % 32; uint32_t warp_id = tid / 32; uint32_t warp_num = blockdim.x / 32; uint32_t offset = blockidx.x * warp_num * hidden_units + warp_id * hidden_units; t val = 0.0f; for (uint32_t i=lane_id; i(d_x, d_y, hidden_units);} 先确定 block_size，这里笔者直接取 128，由于是一个 warp 处理一行数据，所以一个 block 可以处理 warp_num 行数据，总共需要 grid_size 个 block。
核函数内部首先计算当前线程所在的 warp 编号 warp_id 用来定位当前处理元素在哪一行，然后确定线程在 warp 内的编号 lane_id 用来定位该线程具体处理那些元素。由于矩阵宽度 hidden_units 实际肯定还是比 32 大的，所以不可能说一个线程只处理一个元素，因此每个线程会处理多个元素，步长为 32，例如当 hidden_units 为 128 时，lane_id = 0 的线程将处理位置为 0、32、64、96 的四个元素，lane_id = 1 的线程将处理位置为 1、33、65、97 的四个元素，以此类推，这个计算过程是没有并行的。循环计算一轮后，对线程束内每个线程的 val 进行束内规约就可以得到一行元素的规约和。
3 块内规约 块内规约，就是在一个线程块内求规约值，通常块内规约会通过束内规约来实现，以下是 nvidia 提供的一个块内规约设备函数。
template __inline__ __device__t blockreducesum(t val){ static __shared__ t shared[32]; int lane = threadidx.x & 0x1f; int wid = threadidx.x >> 5; val = warpreducesum(val); if(lane == 0) shared[wid] = val; __syncthreads(); val = (threadidx.x > 5 )) ? shared[lane] : (t)0.0f; val = warpreducesum(val); return val;} 规约思路分为两步，首先通过束内规约求出当前线程所在 warp 的规约值，存入 shared 中，然后把 warpsum 赋值给 threadidx.x 小于 32 的线程内的变量 val，这 32 个线程正好也在一个线程束内，然后再执行一次束内规约就得到块内规约值，计算思路非常巧妙。
另外针对块内规约的问题，官方 cub 库其实提供了 api，开发者可以导入头文件 cub/cub.cuh 后直接使用，注意低版本的 cuda 不支持此 api。我们来看下 api 的调用方式。
#include templatestruct sumop { __device__ __forceinline__ t operator()(const t& a, const t& b) const { return a + b; }};template__inline__ __device__ t blockallreduce(t val) { typedef cub::blockreduce blockreduce; __shared__ typename blockreduce::tempstorage temp_storage; __shared__ t result_broadcast; t result = blockreduce(temp_storage).reduce(val, reductionop()); if (threadidx.x == 0) { result_broadcast = result; } __syncthreads(); return result_broadcast;} 除了必要的待规约变量、block_size 以外，还需要传入一个计算函数，笔者给出了示例 sumop。
当矩阵宽度 hidden_units 较大时，通常可以使用一个 block 处理一行数据，笔者给出具体的核函数如下：
template__global__ void matrix2dblockreducesum(const t* inp, t*out, const uint32_t hidden_units) { t val = 0.0f; uint32_t offset = blockidx.x * hidden_units; for (uint32_t i=threadidx.x; i(d_x, d_y, hidden_units);} 同样，block_size 这里笔者直接取 128，由于是一个 block 处理一行数据，总共需要 batch_size 个 block。
由于矩阵宽度 hidden_units 实际肯定还是比 block_size 大的，所以不可能说一个线程只处理一个元素，因此每个线程会处理多个元素，步长为 block_size，例如当 hidden_units 为 512 时，lane_id = 0 的线程将处理位置为 0、128、256、384 的四个元素，lane_id = 1 的线程将处理位置为 1、129、257、385 的四个元素，以此类推，这个计算过程是没有并行的。循环计算一轮后，对 block 内每个线程的 val 进行块内规约就可以得到一行元素的规约和。
4 向量化数据提升访存带宽 使用向量化操作能够提升内存读写的带宽，而 cuda 里也提供了一系列数据类型来支持向量化操作，如 float2、float4，就是将 2 个或 4 个 float 数据作为一个整体。为了增加代码的复用性，笔者这里封装了一个 packed 数据结构，用于对不同的数据类型进行打包。
template struct alignas(sizeof(t) * pack_size) packed{ __device__ packed() { // do nothing } union { t elem[pack_size]; // 这里联合体只有一个成员，为了方便后期扩展 };}; 结构体内有一个 elem 数组变量，整个结构的内存对齐设置为 sizeof(t) * pack_size，说白了其实就是把 pack_size 个 t 类型的数据“捆绑”在一起组成一个新的数据结构，读写内存的时候只需要一次读写就可以读 pack_size 个数据，目的是减小内存读写次数。
那么这个 pack_size 能不能无限大呢？显然不能，cuda 里最大支持 128 bit 的访问粒度，也就是说对于 float 类型（占 4 个字节，32 bit），一次最多读写 4 个，也就是说 float 的 pack_size 最多取到 4，本文笔者的示例代码中数据类型都以 float 为例，pack_size 取 4。
4.1 pack 后的束内规约示例代码 将 matrix2dwarpreducesum 改写为 pack 版的核函数也很简单，计算思路都是一致的，只不过原来一次访问一个元素，现在一次访问一个 pack 的元素，在执行核函数之前笔者加了一个断言，保证 hidden_units 能够被 pack_size 整除，具体代码如下。
template __global__ void matrix2dwarpreducesumpack(const t* d_x, t* d_y, const uint32_t hidden_units, const uint32_t num_packs) { const uint32_t warp_id = threadidx.x / 32; const uint32_t lane_id = threadidx.x & 0x1f; const uint32_t warp_num = blockdim.x / 32; const uint32_t offset = blockidx.x * warp_num * hidden_units + warp_id * hidden_units; const packed* buf = reinterpret_cast(d_x + offset); packed pack; t val = 0.0f; for (uint32_t pack_id=lane_id; pack_id(d_x, d_y, hidden_units, num_packs);} 核函数内部就一句核心代码，将 const t* 指针转换成 const packed*。
const packed* buf = reinterpret_cast(d_x + offset); 然后用 pack_id 索引一次取一个 pack 的数据，注意这里对 pack 索引的时候不要写错了。跟前面一样，相邻的线程处理相邻的 pack 数据，这是为了全局内存的合并访问。加法计算次数还是那么多次，因为 packed 结构体并不能直接参与计算，还是要用 elem 里面的元素计算，这个核函数也就节省了访存次数而已。
4.2 pack 后的块内规约示例代码 matrix2dblockreducesumpack 核函数的实现就更简单了，直接上代码。
template __global__ void matrix2dblockreducesumpack(const t* d_x, t* d_y, const uint32_t hidden_units, const uint32_t num_packs) { t val = 0.0f; uint32_t offset = blockidx.x * hidden_units; const packed* buf = reinterpret_cast(d_x + offset); packed pack; for (uint32_t pack_id=threadidx.x; pack_id(d_x, d_y, hidden_units, num_packs);} 5 小结 在深度学习算子的开发过程中，规约是一个非常常见的场景，以 softmax 为例就有 reducemax 和 reducesum 的应用，本文给出了两种规约实现方式，可供读者参考使用。实际开发过程中，规约计算一般是隐藏在其他 kernel 中的，并不会奢侈到单独写个规约 kernel，所以要求开发人员领会思路活学活用。
Flyme AR系统更新增加提词器与眼镜音乐续播功能
深圳新闻网：智能传感器产业在光明迎来发展的春天
nubiaZ9拆解 无边框到底是如何做到的
板子上面的这些处理器你都知道吗？
低压断路器的三段式保护整定
束内规约与块内规约问题
毕业真实的版本*「莫道克大学毕业证书」Murdoch原件一模一样证书
360视觉云探路中小微企业数智转型 站在数字经济浪尖 驶向数字化
4.8英寸屏双核 Vega Racer 2本月底推出
互补色和色环
AC312E直升机已通过了中国民航局的型号审查
如何用Python开发OpenHarmony设备程序
苹果打算推出三摄手机来对抗华为P20 Pro？
建筑防雷接地、工作接地、保护接地系统的设计
继6月取得了7.4%的增长后，美国的7月制造业产出增长了3.4%
为什么服务器要在机房的托管下才能更好的运行
惠普Sonata锂离子长效电池超长寿命电池开始销售
冠层快速分析仪的作用是怎样的
微软Surface系列的全新产品介绍
通用重整旗鼓向电动车进发 但却夹杂着英雄迟暮的悲壮
1 应用背景 关于规约的定义，相信能读到这篇文章的读者都不陌生，笔者在早期的文章中也介绍过一些规约方法，基本思想都是折半规约，主要应用于较大元素规模的向量规约，有兴趣的读者可以移步【cuda编程】cuda编程中的并行规约问题。
本文要介绍的规约场景与之前有所不同，主要应用于矩阵规约，也就是说本文假设的输入变量的维度是 2 维的，形状为 [batch_size, hidden_units]，规约之后的输出变量形状为 [batch_size, ]。
接下来，本文将以规约求和为例介绍两种规约方式：束内规约、块内规约。
2 束内规约 束内规约，也就是在一个线程束内对某个变量进行规约。我们知道 cuda 架构下指令是以线程束（相邻的 32 个线程）为基本单元执行的，线程束内也可以通过束内洗牌指令进行通信，所以这提供了一个很好的束内规约思路。下面是 nvidia 提供的基础的一个规约设备函数。
template __inline__ __device__t warpreducesum(t val){ for(int mask = 16; mask > 0; mask >>= 1) val += __shfl_xor_sync(final_mask, val, mask, 32); return val;} 这个设备函数可以求出当前线程所在线程束的指定变量的规约和，原理涉及洗牌指令的计算逻辑，不再赘述。
当矩阵宽度 hidden_units 较小时，通常可以使用一个 warp 处理一行数据，一个 block 内可以处理多行数据，笔者给出具体的核函数如下：
// 一个 warp 处理一行数据template__global__ void matrix2dwarpreducesum(const t* inp, t*out, const uint32_t hidden_units) { uint32_t tid = threadidx.x; uint32_t lane_id = tid % 32; uint32_t warp_id = tid / 32; uint32_t warp_num = blockdim.x / 32; uint32_t offset = blockidx.x * warp_num * hidden_units + warp_id * hidden_units; t val = 0.0f; for (uint32_t i=lane_id; i
核函数内部首先计算当前线程所在的 warp 编号 warp_id 用来定位当前处理元素在哪一行，然后确定线程在 warp 内的编号 lane_id 用来定位该线程具体处理那些元素。由于矩阵宽度 hidden_units 实际肯定还是比 32 大的，所以不可能说一个线程只处理一个元素，因此每个线程会处理多个元素，步长为 32，例如当 hidden_units 为 128 时，lane_id = 0 的线程将处理位置为 0、32、64、96 的四个元素，lane_id = 1 的线程将处理位置为 1、33、65、97 的四个元素，以此类推，这个计算过程是没有并行的。循环计算一轮后，对线程束内每个线程的 val 进行束内规约就可以得到一行元素的规约和。
3 块内规约 块内规约，就是在一个线程块内求规约值，通常块内规约会通过束内规约来实现，以下是 nvidia 提供的一个块内规约设备函数。
template __inline__ __device__t blockreducesum(t val){ static __shared__ t shared[32]; int lane = threadidx.x & 0x1f; int wid = threadidx.x >> 5; val = warpreducesum(val); if(lane == 0) shared[wid] = val; __syncthreads(); val = (threadidx.x > 5 )) ? shared[lane] : (t)0.0f; val = warpreducesum(val); return val;} 规约思路分为两步，首先通过束内规约求出当前线程所在 warp 的规约值，存入 shared 中，然后把 warpsum 赋值给 threadidx.x 小于 32 的线程内的变量 val，这 32 个线程正好也在一个线程束内，然后再执行一次束内规约就得到块内规约值，计算思路非常巧妙。
另外针对块内规约的问题，官方 cub 库其实提供了 api，开发者可以导入头文件 cub/cub.cuh 后直接使用，注意低版本的 cuda 不支持此 api。我们来看下 api 的调用方式。
#include templatestruct sumop { __device__ __forceinline__ t operator()(const t& a, const t& b) const { return a + b; }};template__inline__ __device__ t blockallreduce(t val) { typedef cub::blockreduce blockreduce; __shared__ typename blockreduce::tempstorage temp_storage; __shared__ t result_broadcast; t result = blockreduce(temp_storage).reduce(val, reductionop()); if (threadidx.x == 0) { result_broadcast = result; } __syncthreads(); return result_broadcast;} 除了必要的待规约变量、block_size 以外，还需要传入一个计算函数，笔者给出了示例 sumop。
当矩阵宽度 hidden_units 较大时，通常可以使用一个 block 处理一行数据，笔者给出具体的核函数如下：
template__global__ void matrix2dblockreducesum(const t* inp, t*out, const uint32_t hidden_units) { t val = 0.0f; uint32_t offset = blockidx.x * hidden_units; for (uint32_t i=threadidx.x; i
由于矩阵宽度 hidden_units 实际肯定还是比 block_size 大的，所以不可能说一个线程只处理一个元素，因此每个线程会处理多个元素，步长为 block_size，例如当 hidden_units 为 512 时，lane_id = 0 的线程将处理位置为 0、128、256、384 的四个元素，lane_id = 1 的线程将处理位置为 1、129、257、385 的四个元素，以此类推，这个计算过程是没有并行的。循环计算一轮后，对 block 内每个线程的 val 进行块内规约就可以得到一行元素的规约和。
4 向量化数据提升访存带宽 使用向量化操作能够提升内存读写的带宽，而 cuda 里也提供了一系列数据类型来支持向量化操作，如 float2、float4，就是将 2 个或 4 个 float 数据作为一个整体。为了增加代码的复用性，笔者这里封装了一个 packed 数据结构，用于对不同的数据类型进行打包。
template struct alignas(sizeof(t) * pack_size) packed{ __device__ packed() { // do nothing } union { t elem[pack_size]; // 这里联合体只有一个成员，为了方便后期扩展 };}; 结构体内有一个 elem 数组变量，整个结构的内存对齐设置为 sizeof(t) * pack_size，说白了其实就是把 pack_size 个 t 类型的数据“捆绑”在一起组成一个新的数据结构，读写内存的时候只需要一次读写就可以读 pack_size 个数据，目的是减小内存读写次数。
那么这个 pack_size 能不能无限大呢？显然不能，cuda 里最大支持 128 bit 的访问粒度，也就是说对于 float 类型（占 4 个字节，32 bit），一次最多读写 4 个，也就是说 float 的 pack_size 最多取到 4，本文笔者的示例代码中数据类型都以 float 为例，pack_size 取 4。
4.1 pack 后的束内规约示例代码 将 matrix2dwarpreducesum 改写为 pack 版的核函数也很简单，计算思路都是一致的，只不过原来一次访问一个元素，现在一次访问一个 pack 的元素，在执行核函数之前笔者加了一个断言，保证 hidden_units 能够被 pack_size 整除，具体代码如下。
template __global__ void matrix2dwarpreducesumpack(const t* d_x, t* d_y, const uint32_t hidden_units, const uint32_t num_packs) { const uint32_t warp_id = threadidx.x / 32; const uint32_t lane_id = threadidx.x & 0x1f; const uint32_t warp_num = blockdim.x / 32; const uint32_t offset = blockidx.x * warp_num * hidden_units + warp_id * hidden_units; const packed* buf = reinterpret_cast(d_x + offset); packed pack; t val = 0.0f; for (uint32_t pack_id=lane_id; pack_id
const packed* buf = reinterpret_cast(d_x + offset); 然后用 pack_id 索引一次取一个 pack 的数据，注意这里对 pack 索引的时候不要写错了。跟前面一样，相邻的线程处理相邻的 pack 数据，这是为了全局内存的合并访问。加法计算次数还是那么多次，因为 packed 结构体并不能直接参与计算，还是要用 elem 里面的元素计算，这个核函数也就节省了访存次数而已。
4.2 pack 后的块内规约示例代码 matrix2dblockreducesumpack 核函数的实现就更简单了，直接上代码。
template __global__ void matrix2dblockreducesumpack(const t* d_x, t* d_y, const uint32_t hidden_units, const uint32_t num_packs) { t val = 0.0f; uint32_t offset = blockidx.x * hidden_units; const packed* buf = reinterpret_cast(d_x + offset); packed pack; for (uint32_t pack_id=threadidx.x; pack_id
Flyme AR系统更新增加提词器与眼镜音乐续播功能
深圳新闻网：智能传感器产业在光明迎来发展的春天
nubiaZ9拆解 无边框到底是如何做到的
板子上面的这些处理器你都知道吗？
低压断路器的三段式保护整定
束内规约与块内规约问题
毕业真实的版本*「莫道克大学毕业证书」Murdoch原件一模一样证书
360视觉云探路中小微企业数智转型 站在数字经济浪尖 驶向数字化
4.8英寸屏双核 Vega Racer 2本月底推出
互补色和色环
AC312E直升机已通过了中国民航局的型号审查
如何用Python开发OpenHarmony设备程序
苹果打算推出三摄手机来对抗华为P20 Pro？
建筑防雷接地、工作接地、保护接地系统的设计
继6月取得了7.4%的增长后，美国的7月制造业产出增长了3.4%
为什么服务器要在机房的托管下才能更好的运行
惠普Sonata锂离子长效电池超长寿命电池开始销售
冠层快速分析仪的作用是怎样的
微软Surface系列的全新产品介绍
通用重整旗鼓向电动车进发 但却夹杂着英雄迟暮的悲壮