当前位置：首页 > news >正文

《CUDA编程》10.线程束的基本函数与协作组

news 2024/10/27 22:07:42

本章将对线程束及相关的 CUDA 语法和编程模型进行更加深入、系统的介绍，涉及的大部分 CUDA 编程功能由 CUDA 9 引入，而且有些功能仅受到帕斯卡或更高架构的支持。

1单指令-多线程执行模式

1.1 在伏特架构之前

在伏特架构之前，一个线程束中的线程拥有同一个程序计时器，但各自有不同的寄存器状态，从而可以根据程序的逻辑判断选择不同的分支，但是在执行的时候是按照顺序执行。即，在同一时刻，一个线程束中的线程只能执行同一个指令，或者都不执行，这就是单指令-多线程（SIMT）的执行模式

当一个线程束中的线程顺序地执行判断语句中的不同分支时，我们称发生了分支发散（branch divergence），例如，假设核函数中有如下判断语句：

if (condition)
{A;
}
else{B;
}

**首先，满足 condition 的线程会执行语句 A，其他的线程将闲置。然后，不满足 condition 的线程会执行语句 B，其他的线程将闲置。**所以如果语句A的指令和语句B的指令数量差不多时，每次执行都有一半的线程闲置，执行效率也比没有分支的情况低了一倍。

注意： $分支发散$ 是针对同一个线程束内部的线程，如果是不同的线程束执行不同的分支，则不属于，观察下面的代码：

int warp_id=threadIdx.x / 32;//32个线程是一个线程束
switch (warp_id) { case 0 : S0; break; case 1 : S1; break; case 2 : S2; break; case 3 : S3; break; 
}

这个switch语句根据warp_id的值来执行不同的代码块：

如果warp_id是0，执行S0
如果warp_id是1，执行S1
如果warp_id是2，执行S2
如果warp_id是3，执行S3

通过这种方式，可以对不同的warp分配不同的任务。由于warp中的线程是同步执行的，这种设计可以让不同的warp高效并行地执行不同的代码块。

再观察下面的代码：

int lane_id = threadIdx.x % 32; // 通过取余，可以获取该线程在当前线程束的位置
switch (lane_id) 
{ case 0 : S0; break; case 1 : S1; break;  case 2 : S2; break; 
}

如果lane_id是0，执行S0
如果lane_id是1，执行S1
如果lane_id是2，执行S2

也就是说每次只有一个线程在执行，其他31个都闲置，将导致严重的分支发散。

所以，在编写核函数时要尽量避免分支发散，但有时又无法避免，比如我们每次在核函数中运行计算算法前，会有一个if(n<N){do something;}的防止越界代码，虽然会影响性能，但是为了程序的正确性，不得不牺牲这点性能。

1.2 从伏特架构开始

从伏特架构开始，引入了独立线程调度（independent thread scheduling）机制。每个线程有自己的程序计数器。这使得伏特架构有了一些以前的架构所没有的新的线程束内同步与通信的模式，从而提高了编程的灵活性。

要实现独立线程调度机制，一个代价是增加了寄存器负担：单个线程的程序计数器一般需要使用两个寄存器。也就是说，伏特架构的独立线程调度机制使得 SM 中每个线程可利用的寄存器少了两个。另外，独立线程调度机制使得假设了线程束同步（warp synchronous）的代码变得不再安全。

例如，在数组归约的例子中，当线程号小于 32 时，可以省去线程块同步函数 __syncthreads();，但是在伏特架构中便不再安全，所以我们将介绍一种比线程块同步函数更细粒的同步函数，线程束内同步函数__syncwarp

2 线程束内的线程同步函数

当所涉及的线程都在一个线程束内时，可以将线程块同步函数 __syncthreads 换成一个更加廉价的线程束同步函数 __syncwarp。我们将它简称为束内同步函数。该函数的原型为:

void __syncwarp(unsigned mask = 0xffffffff);

该函数有一个可选的参数。该参数是一个代表掩码的无符号整型数，默认值的全部 32 个二进制位都为 1，代表线程束中的所有线程都参与同步,例如，掩码 0xfffffffe 代表排除第 0 号线程。用该束内同步函数。下面是将之前的归约求和核函数进行改写：

__global__ void reduce_syncwarp(const float* d_x, float* d_y, const int N) {const int tid = threadIdx.x;const int bid = blockIdx.x;const int n = bid * blockDim.x + tid;extern __shared__ float s_y[];s_y[tid] = (n<N)? d_x[n]:0.0f;__syncthreads();for (int offset = blockDim.X >> 1; offset >= 32; offset >>= 1) {if (tid < offset) {s_y[tid] += s_y[tid + offset];}__syncthreads();}for (int offset = 16; offset > 0; offset >>= 1) {if (tid < offset) {s_y[tid] += s_y[tid + offset];}__syncwarp();}if (tid == 0) {atomicAdd(d_y, y);}
}

当offset<=16时，会使用__syncwarp();进行线程内同步，避免不必要的同步，不过即使不使用，也能自动同步，不过会增加不必要的同步开销，但从伏特架构开始，建议时显式的使用，不然可能变得不安全。

当使用共享内存和__syncwarp();时，容易犯一个错误。例如，我们取消了当offset<=16时，取消if (tid < offset){}的判断，假设offset=16，就会出现下列两个计算

s_y[0] += s_y[16]; 
s_y[16] += s_y[32];

即又要读取s_y[16]，又要写入s_y[16]，产生了读写竞争，所以需要一个判断语句来排除该风险

3 更多线程束内的基本函数

这里只介绍线程束表决函数（Warp Vote Functions）和线程束洗牌函数（Warp Shuffle Functions），线程束表决函数和线程束洗牌函数是CUDA中用于warp级别（通常为32个线程）的高级操作，专门用于warp内线程的高效通信和数据处理。它们可以减少共享内存的使用，提高warp内部数据传输和逻辑判断的效率，进而优化性能。

①表决函数主要用于在warp内的多个线程之间进行逻辑判断，以判断warp内特定条件是否满足，线程束表决函数的原型如下：

unsigned __ballot_sync(unsigned mask, int predicate);
int __all_sync(unsigned mask, int predicate);
int __any_sync(unsigned mask, int predicate);

②线程束洗牌函数用于warp内线程之间直接交换数据，无需借助共享内存。这种warp内数据交换可以显著降低内存访问延迟，并且使用起来更为高效，线程束洗牌函数的原型如下：

T __shfl_sync(unsigned mask, T v, int srcLane, int w = warpSize);
T __shfl_up_sync(unsigned mask, T v, unsigned d, int w = warpSize);
T __shfl_down_sync(unsigned mask, T v, unsigned d, int w = warpSize);
T __shfl_xor_sync(unsigned mask, T v, int laneMask, int w =
warpSize);

其中，类型 T 可以为整型（int）、长整型（long）、长长整型（long long）、无符号整型（unsigned）、无符号长整型（unsigned long）、无符号长长整型（unsigned long long）、单精度浮点型（float）及双精度浮点型（double）。

这些函数的功能如下：

__ballot_sync(mask, predicate)：该函数返回一个无符号整数。如果线程束内第 n 个线程参与计算且 predicate 值非零，则将所返回无符号整数的第 n 个二进制位取为1，否则取为 0。这里，参与的线程对应于 mask 中取1的比特位。该函数的功能相当于从一个旧的掩码出发，产生一个新的掩码。
__all_sync(mask, predicate)：线程束内所有参与线程的 predicate 值都不为零才返回 1，否则返回 0。这里参与的线程对应于 mask 中取 1 的比特位。该函数实现了一个"规约-广播"（reduction-and-broadcast）式计算。该函数类似于这样一种选举操作：当所有参选人都同意时才通过。
__any_sync(mask, predicate,)：线程束内所有参与线程的 predicate 值有一个不为零就返回 1，否则返回 0。这里参与的线程对应于 mask 中取 1 的比特位。该函数也实现了一个"规约-广播"式计算。该函数类似于这样一种选举操作：只要有一个参选人同意就通过。
__shfl_sync(mask, v, srcLane, w)：参与线程返回标号为 srcLane 的线程中变量 v 的值。这是一种广播式数据交换：将一个线程中的数据广播到所有（包括自己）线程。
__shfl_up_sync(mask, v, d, w)：标号为 t 的参与线程返回标号为 t - d 的线程中变量 v 的值。标号满足 t - d < 0 的线程返回原来的 v。例如：当 w = 8， d = 2 时，该函数将第 0-5 号线程中变量 v 的值传送到第 2-7 号线程，而第 0-1 号线程返回它们原来的 v。形象地说，这是一种将数据向上平移的操作。
__shfl_down_sync(mask, v, d, w)：标号为 t 的参与线程返回标号为 t + d 的线程中变量 v 的值。标号满足 t + d >= w 的线程返回原来的 v。例如：当 w = 8，d = 2 时，该函数将第 2-7 号线程中变量 v 的值传送到第 0-5 号线程，而第 6-7 号线程返回它们原来的 v。形象地说，这是一种将数据向下平移的操作。
__shfl_xor_sync(mask, v, laneMask, w)：标号为 t 的参与线程返回标号为 $t ^ {laneMask}$ 的线程中变量 v 的值。这里， $t ^ {laneMask}$ 表示两个整数按位异或运算的结果。
例如，当 w = 8， laneMask = 2 时，第 0-7 号线程的按位异或运算 $t ^ {laneMask}$ 分别如下：
0 ^ 2 = 0000 ^ 0010 = 0010 = 2
1 ^ 2 = 0001 ^ 0010 = 0011 = 3
2 ^ 2 = 0010 ^ 0010 = 0000 = 0
3 ^ 2 = 0011 ^ 0010 = 0001 = 1
4 ^ 2 = 0100 ^ 0010 = 0110 = 6
5 ^ 2 = 0101 ^ 0010 = 0111 = 7
6 ^ 2 = 0110 ^ 0010 = 0100 = 4
7 ^ 2 = 0111 ^ 0010 = 0101 = 5

我们不需要在任何地方明显地使用同步函数，如__syncwarp。这是因为，这里所有的线程束内的基本函数（都以 _sync 结尾）都具有隐式的同步功能。从伏特架构开始，在使用这些函数时必须使用由 CUDA 9 引入的新版本，不要再使用原来的没有 _sync 的版本。

3.1 利用线程束洗牌函数进行归约计算

__global__ void reduce_syncwarp(const float* d_x, float* d_y, const int N) {const int tid = threadIdx.x;const int bid = blockIdx.x;const int n = bid * blockDim.x + tid;extern __shared__ float s_y[];s_y[tid] = (n<N)? d_x[n]:0.0f;__syncthreads();for (int offset = blockDim.X >> 1; offset >= 32; offset >>= 1) {if (tid < offset) {s_y[tid] += s_y[tid + offset];}__syncthreads();}for (int offset = 16; offset > 0; offset >>= 1) {y += __shfl_down_sync(FULL_MASK, y, offset);}if (tid == 0) {atomicAdd(d_y, y);}
}

假设 blockDim.x 为 64，那么 offset 的变化过程如下：

共享内存归约阶段（offset >= 32）：

offset = 32: tid < 32 的线程将 s_y[tid] += s_y[tid + 32]
offset = 16: tid < 16 的线程将 s_y[tid] += s_y[tid + 16]
offset = 8: tid < 8 的线程将 s_y[tid] += s_y[tid + 8]
offset = 4: tid < 4 的线程将 s_y[tid] += s_y[tid + 4]
offset = 2: tid < 2 的线程将 s_y[tid] += s_y[tid + 2]
offset = 1: tid < 1 的线程将 s_y[tid] += s_y[tid + 1]

warp shuffle 归约阶段（offset <= 16）：

offset = 16: 使用 __shfl_down_sync，每个线程从相距16的邻居读取数据并相加。
offset = 8: 每个线程从相距8的邻居读取数据并相加。
offset = 4: 每个线程从相距4的邻居读取数据并相加。
offset = 2: 每个线程从相距2的邻居读取数据并相加。
offset = 1: 每个线程从相距1的邻居读取数据并相加。

相比之前的版本有两处不同：
第一，在进行线程束内的循环之前，这里将共享内存中的数据复制到了寄存器。在线程束内使用洗牌函数进行规约时，不再需要明显地使用共享内存。

第二，用语句y += __shfl_down_sync(FULL_MASK, y, offset);替换了语句块if (tid < offset){}，去掉了同步函数，也去掉了对线程号的限制，因为洗牌函数能够自动处理同步与读-写竞争问题。

4 协作组

协作组（Cooperative Groups）是一种抽象，用于组织和管理线程之间的协作和同步。协作组通过提供灵活的线程分组和同步方式，使线程能够在更小的范围内协同工作，而不仅限于整个线程块或warp级别的协作。

协作组主要通过 cooperative_groups 命名空间中的API提供各种线程分组和操作方式。它允许开发者根据需求，将线程分为不同的协作组，每个协作组中的线程可以同步执行任务或共享数据。

在CUDA中，协作组的基本单位包括：

线程块级协作组（Block-wide Group）：包括整个线程块内的所有线程。
warp级协作组（Warp-wide Group）：包括一个warp（通常是32个线程）内的所有线程。
线程子组（Thread Subgroups）：可以是任意数量的线程组合，灵活性更高。
线程组内的单一线程（Single Thread Group）：可以是一个单独的线程，常用于需要单线程执行的操作。

使用协作组的功能时需要在相关源文件包含如下头文件：
#include <cooperative_groups.h>

除此以外，还要注意所有与协作组相关的数据类型和函数都定义在名字空间（namespace）cooperative_groups 下。可以用如下语句导入该名字空间中的内容：using namespace cooperative_groups;也可以给该名字空间起一个较短的别名：namespace cg = cooperative_groups;

4.1 线程块级别的协作组

协作组编程模型中最基本的类型是线程组 thread_group，该类型有如下成员：

void sync(); 该函数能同步组内所有线程。
unsigned size(); 该函数返回组内总的线程数目，即组的大小。
unsigned thread_rank(); 该函数返回当前调用该函数的线程在组内的标号（从 0开始计数）。
bool is_valid(); 该函数返回一个逻辑值，如果定义的组违反了任何 CUDA 的限制，则为假，否则为真。

线程组类型有一个称为线程块 thread_block 的导出类型，可以用thread_block g = this_thread_block();定义并初始化一个 thread_block 对象，在该类型中提供了两个额外的函数：

dim3 group_index(); 该函数返回当前调用该函数的线程的线程块指标，等价于 blockIdx。
dim3 thread_index(); 该函数返回当前调用该函数的线程的线程指标，等价于 threadIdx。

其中，this_thread_block() 相当于一个线程块类型的常量。
这样定义的 g 就代表我们已经非常熟悉的线程块，只不过这里把它包装成了一个类型。g.sync() 完全等价于 __syncthreads()，g.group_index() 完全等价于 CUDA 中的内建变量 blockIdx，g.thread_index() 完全等价于 CUDA 中的内建变量 threadIdx。

可以用函数 tiled_partition 将一个线程块划分为若干片（tile），每一片构成一个新的线程组。
目前仅仅可以将片的大小设置为 2 的正整数次方且不大于 32，也就是 2、4、8、16 和 32（和线程束洗牌函数的最后一个参数类似）。例如，如下语句通过函数 tiled_partition 将一个线程块分割为我们熟知的线程束：thread_group g32 = tiled_partition(this_thread_block(), 32);

我们还可以将该线程组分割为更细的线程组。如下语句将每个线程束再分割为包含4个线程的线程组：thread_group g4 = tiled_partition(g32, 4);

当这种线程组的大小在编译期间就已知时，可以用如下模板化的版本（可能更加高效）进行定义：

thread_block_tile<32> g32 = tiled_partition<32>(this_thread_block());
thread_block_tile<4> g4 = tiled_partition<4>(this_thread_block());

这样定义的线程组一般称为线程块片（thread block tile）。线程块片还额外地定义了如下函数（类似于线程束内的基本函数）：

unsigned __ballot_sync(int predicate);
int __all_sync(int predicate);
int __any_sync(int predicate);
T __shfl_sync(T v, int srcLane);
T __shfl_up_sync(T v, unsigned d);
T __shfl_down_sync(T v, unsigned d);
T __shfl_xor_sync(T v, int laneMask);

线程块片的函数有两点不同：
第一，线程块片的函数少了第一个代表掩码的参数，因为线程组内的所有线程都必须参与相关函数的运算。

第二，线程块片的洗牌函数（上述函数中的后 4 个）少了最后一个代表宽度的参数，因为该宽度就是线程块片的大小，即定义线程块片的模板参数。

4.2 利用协作组进行归约计算

__global__ void reduce_syncwarp(const float* d_x, float* d_y, const int N) {const int tid = threadIdx.x;const int bid = blockIdx.x;const int n = bid * blockDim.x + tid;extern __shared__ float s_y[];s_y[tid] = (n<N)? d_x[n]:0.0f;__syncthreads();for (int offset = blockDim.X >> 1; offset >= 32; offset >>= 1) {if (tid < offset) {s_y[tid] += s_y[tid + offset];}__syncthreads();}float y = s_y[tid];// 创建一个warp级协作组，大小为32thread_block_tile<32> g = tile_partition<32>(this_thread_block());for (int i = g.size() >> 1; i > 0; i >>= 1) {y += g.shfl_down(y, i);// 从偏移为i的线程中读取数据并相加}if (tid == 0) {atomicAdd(d_y, y);}
}

优化到这一步，归约计算模块语句相比于CPU已经快了40倍。

5 数组归约程序的进一步优化

5.1 提高线程利用率

在 $4$ 中的归约核函数中，其实线程的利用率并不高。

因为我们使用大小为 128 的线程块：

当 offset 等于 64 时，只用了 1/2 的线程进行计算，其余线程闲置。
当 offset 等于 32 时，只用了 1/4 的线程进行计算，其余线程闲置。
当 offset 等于 1 时，只用了 1/128 的线程进行计算，其余线程闲置。

相比之下，在归约之前，将全局内存中的数据复制到共享内存的操作（见第 7 行）对线程的利用率是 100% 的。

据此得到一个想法：如果能够提高归约之前所做计算的比例，那应该可以从整体上提升对线程的利用率（在这一过程中，除了简单的搬运数据外，我们能够进行一些额外的计算，那么这段时间内，线程的工作负载就会增加，从而提升线程的利用率）。

__global__ void reduce_cp(const float* d_x, float* d_y, const int N) {const int tid = threadIdx.x;const int bid = blockIdx.x;extern __shared__ float s_y[];float y = 0.0;const int stride = blockDim.x*gridDim.x;//步长，确保每个线程处理不同的数据元素。for (int n = bid * blockDim.x + tid; n < N; n += stride) {y += d_x[n];//循环遍历输入数组 d_x，每个线程负责累加自己负责的数据片段。}s_y[tid] = y;//将每个线程计算的局部和存入共享内存，并进行同步__syncthreads();for (int offset = blockDim.X >> 1; offset >= 32; offset >>= 1) {if (tid < offset) {s_y[tid] += s_y[tid + offset];}__syncthreads();}y = s_y[tid];thread_block_tile<32> g = tile_partition<32>(this_thread_block());for (int i = g.size() >> 1; i > 0; i >>= 1) {y += g.shfl_down(y, i);}if (tid == 0) {d_y[bid] = y;}
}

5.2 避免反复分配与释放设备内存

下面给出一个调用该核函数的包装函数，它将返回最终的结果。这里，我们将GRID_SIZE取为10240，将BLOCK_SIZE取为128：

float reduce(const float* d_x) {const int ymem = sizeof(float) * GRID_SIZE;const int smem = sizeof(float) * BLOCK_SIZE;float h_y[1] = { 0 };float* d_y;CHECK(cudaMalloc(&d_y, ymem));reduce_cp << <GRID_SIZE, BLOCK_SIZE, smem >> > (d_x, d_y, N);reduce_cp << <1, 1024, sizeof(real) * 1024 >> > (d_y, d_y, GRID_SIZE);CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDeviceToHost));CHECK(cudaFree(d_y));return h_y[0];
}

这一版代码相较于上一版最优的代码，又提升了40%

在上面的包装函数 reduce 中，我们需要为数组 d_y 分配与释放设备内存。实际上，设备内存的分配与释放是比较耗时的。一种优化方案是使用静态全局内存代替这里的动态全局内存，因为静态内存是编译期间就会分配好的，不会在运行程序时反复地分配，故比动态内存分配高效很多。

可以用如下语句在函数外部定义我们需要的静态全局内存变量：device real static_y[GRID_SIZE];

我们可以直接在核函数中使用该变量，但这需要改变核函数代码。如果不想改变核函数代码，可以利用运行时 API 函数 cudaGetSymbolAddress 获得一个指向该静态全局内存的指针，供核函数使用，该函数的原型如下：

cudaError_t cudaGetSymbolAddress(void **devPtr, const void *symbol);

这里的 symbol 参数可以是静态全局内存（用 device 定义）或者常量内存（用 constant 定义）的变量名。通过函数 cudaGetSymbolAddress 获得的设备指针可以像其他设备指针一样使用。修改后函数如下：

float reduce(const float* d_x) {float* d_y;CHECK(cudaGetSymbolAddress((void**)&d_y, static_y));const int smem = sizeof(real) * BLOCK_SIZE;reduce_cp << <GRID_SIZE, BLOCK_SIZE, smem >> > (d_x, d_y, N);reduce_cp << <1, 1024, sizeof(real) * 1024 >> > (d_y, d_y, GRID_SIZE);real h_y[1] = { 0 };CHECK(cudaMemcpy(h_y, d_y, sizeof(real), cudaMemcpyDeviceToHost));return h_y[0]; 
}