线程束基本函数与协作组

参考书目： CUDA编程基础与实践

作者: 樊哲勇

章节：第10章线程束基本函数和协作组

SIMT单指令-多线程执行模式¶

一个线程束中的多个线程执行相同的指令

分支发散 branch divergence¶

一个线程束中的线程顺序执行判断语句的不同分支时会发生分支发散,只针对一个线程束中的线程

下面这段代码

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

满足condition的线程会先执行A,其他线程会被闲置,执行完之后才轮到不满足condition的线程执行B,其他线程闲置,所以要尽量避免分支发散,但这不可避免,比如判断tid<N,所以要先保证程序正确

独立线程调度 independent thread scheduling¶

从伏特架构开始,引入了独立线程调度机制,每个线程有自己的程序计数器,使其线程间可以同步与通信,代价是增加了寄存器负担,单个线程的程序计数器需要两个寄存器,也就是说,这种独立线程调度机制使得SM种的每个线程可利用的寄存器少了两个,另外,独立线程调度机制没有假设线程束同步,也即是说你需要自己指定线程束同步,也就是比线程块同步函数__syncthreads()颗粒度更低的线程束同步函数__syncwarp,假如不想再伏特架构的GPU上使用独立线程调度机制,就需要在编译时指定虚拟架构低于伏特架构的计算能力,也就是

-arch=compute_60 -code=sm_70

在归约问题中,若所有线程位于一个线程束中,我们就可以使用__syncwarp()代替__syncthreads(),其原型为

void __syncwarp(unsigned mask = 0xfffffff)

可选参数为掩码,其32个二进制数对应着线程束中32个线程,0/1为其对应的线程是否参与同步,默认全部参与, 所以归约核函数更改为

void __global__ reduce_syncwarp(const real *d_x, real *d_y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int n = bid * blockDim.x + tid;
    extern __shared__ real s_y[];
    s_y[tid] = (n < N) ? d_x[n] : 0.0;
    __syncthreads();//否则容易s_y还没初始化好就开始累加了

    for (int offset = blockDim.x >> 1; offset >= 32; offset >>= 1)//如果计算的长度一半大于等于32
    {
        if (tid < offset)//只计算前面一半的线程
        {
            s_y[tid] += s_y[tid + offset];
        }
        __syncthreads();//否则容易到了下个循环,上个循环的线程还没计算完
    }

    for (int offset = 16; offset > 0; offset >>= 1)//计算长度的一半小于32,也就是16了,就可以看作一个线程束
    {
        if (tid < offset)
        {
            s_y[tid] += s_y[tid + offset];
        }
        __syncwarp();//同步此线程束,而不用同步整个线程块
    }

    if (tid == 0)
    {
        atomicAdd(d_y, s_y[0]);
    }
}

分析下面的情况

for (int offset = 16; offset > 0; offset >>= 1)/
{
    s_y[tid] += s_y[tid + offset];
    __syncwarp();//同步此线程束
}

if (tid < offset),制定了一些线程运行,其他线程限制,如果把此限制解除,就是所有线程一起运行,看起来减少了代码量,更加简洁,但是会引发读写竞争, 也就是比如这两个线程s_y[0] += s_y[16]和s_y[16] += s_y[32],同时对s_y[16],从而可能导致结果的错误,所以应该为

real v = 0;
for (int offset = 16; offset > 0; offset >>= 1)/
{
    v += s_y[tid + offset];
    __syncwarp();//同步此线程束(读取操作)
    s_y[tid] += s_y[tid + offset];
    __syncwarp();//同步此线程束(写操作)
}

线程束内基本函数¶

线程表决函数warp vote function¶

unsigned __ballot_sync(unsigned mask, int predicate);

由旧掩码生成新掩码, 在参与计算(mask)的线程中,返回那些predicate为真(非零)的线程(新的mask)

int __all_sync(unsigned mask, int predicate);

在参与计算的线程中,predicate全为真,就返回1,否则返回0,也就是所有参选人都同意,才能返回1

int __any_sync(unsigned mask, int predicate);

在参与计算的线程中,predicate只要有一个为真,就返回1,否则返回0,也就是所有参选人至少一人同意,就返回1

线程洗牌函数¶

这里提出了一个束内指标,通俗可以理解为在线程束内继续分块,由参数w控制,其只能取2,4,6,8,16,32,默认为32,也就是一个束内块中由w个线程

int lane_id = threadIdx.x%w;

由于w为2的幂次,可以用更高效的位运算(按位与)代替取模

int lane_id = threadIdx.x&(w-1);

假设线程块大小为16,假设w=8

线程id:0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
束内id:0 1 2 3 4 5 6 7 0 1 2  3  4  5  6  7

T __shfl_sync(unsigned mask, T v, int srcLane, int w = warpSize)

返回束内块中束内指标为srcLane的v,这是一种广播机制

T __shfl_up_sync(unsigned mask, T v, unsigned d, int w = warpSize)

束内块中,束内指标为t的返回t-d线程中的变量v,t-d<0时,就返回其本身,也就是向上平移操作

T __shfl_down_sync(unsigned mask, T v, unsigned d, int w = warpSize)

束内块中,束内指标为t的返回t+d线程中的变量v,t+d>=w时,返回本身,就也就是向下平移操作

T __shfl_xor_sync(unsigned mask, T v, int laneMask, int w = warpSize)

束内块中,束内指标为t的线程返回指标为t^laneMaske(t按位异或laneMask)线程的变量v,通俗将就是将threadIdx.x与threadIdx.x+laneMask线程之间互换数据.

以上函数均以_sync结尾,就代表了其都具有线程束同步功能,故不需要调用__sync_warp

归约计算¶

使用线程洗牌函数对其进行归约

void __global__ reduce_shfl(const real *d_x, real *d_y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int n = bid * blockDim.x + tid;
    extern __shared__ real s_y[];
    s_y[tid] = (n < N) ? d_x[n] : 0.0;
    __syncthreads();

    for (int offset = blockDim.x >> 1; offset >= 32; offset >>= 1)
    {
        if (tid < offset)
        {
            s_y[tid] += s_y[tid + offset];
        }
        __syncthreads();
    }

    real y = s_y[tid];//下边省略的判断，所以要先把s_y[tid]读出来

    for (int offset = 16; offset > 0; offset >>= 1)
    {
        y += __shfl_down_sync(FULL_MASK, y, offset);
    }

    if (tid == 0)
    {
        atomicAdd(d_y, y);
    }
}

__shfl_down_sync(FULL_MASK, y, offset)表示了所有线程向上平移offset的线程对应的y值,也就等价于y[tid + offset]

同样可以使用__shfl_xor_sync()完成同样的效果,__shfl_xor_sync(FULL_MASK, y, offset)返回值为该线程束内指标+offset的值,也就等价于y[tid + offset] 下面是例子offset=4时

y:               0,1,2,3,4,5,6,7
__shfl_down_sync:4,5,6,7,4,5,6,7
__shfl_xor_sync: 4,5,6,7,0,1,2,3

这样做相比于之前的做法,不用显示的使用共享内存,而是__shfl_down_sync()将其赋值进了寄存器,速度也就更高效,同时去掉了同步函数,自动处理读写竞争问题

协作组¶

在有些并行算法中,需要若干线程间的协作,协作组可以为其提供灵活的协作方式,涉及到线程块内部协作,线程块之间协作,多设备之间协作,这里只讲解线程块内部写作

#include<cooperative_groups.h>//头文件
using namespace cooperative_groups;

线程组thread_group的成员

void sync() 同步组内线程
unsigned size()返回组内线程数
unsigned thread_rank()返回当前线程在组内的标号
bool is_valid()是否违反CUDA限制

其中有成员变量thread_block,包含两个函数

dim3 group_index() 返回当前线程块索引,相当于blockIdx.x
dim3 thread_index() 返回当前线程索引,相当于threadIdx.x

thread_block g = this_thread_block();
g.sync();//等价于__syncthreads()
g.group_index();//等价于blockIdx.x

tiled_partition() 可以将一个线程块划分为若干片,每一片构成一个新的线程组,可以将片大小设置为2的幂次方且不超过32,与洗牌函数的w类似

我们可以通过下面语句将一个线程块分割为线程束

thread_group g32 = tiled_partition(this_thread_block(),32);
//将线程组更细分,一个组有4个线程
thread_group g4 = tiled_partition(g32, 4);
//在编译时就已知线程组大小可用模板化的版本
thread_block_tile<32> g32 = tile_partition<32>(this_thread_block());

这也叫做线程片,有和线程束内基本函数类似的表决和洗牌函数,但是不需要mask和w参数,默认全部参与计算

归约计算¶

利用协作组归约计算

y = s_y[tid];
thread_block_tile<32> g = tiled_partition<32>(this_thread_block());
for (int i = g.size() >> 1; i > 0; i >>= 1)
{
    y += g.shfl_down(y, i);
}
// 类似线程束内函数的实现方式，这时线程束就是32
// for (int offset = 16; offset > 0; offset >>= 1)
// {
//     y += __shfl_down_sync(FULL_MASK, y, offset);
// }

对归约函数进一步优化¶

我们之前的例子可以看出，我们只计算offset以内的线程，每经过一次循环，offset减半，参与计算的线程数也就减半，其余线程闲置，但是全局内存到共享内存线程利用率时100%，可以在归约之间加大计算比例，也就是不之间将全局内存复制到共享内存中，而是对其先进行一部分累加，共享内存存累加后的结果

为了做到这一点，一个线程要处理若干数据，要注意的时，一个线程处理的数据要尽可能间隔尽量远(一个线程块或一个网格)，因为相邻的数据要给相邻的线程计算，这样才能保证全局内存的合并访问

我们此时的例子并不满足GRID_SIZE = (N+BLOCK_SIZE-1)/BLOCK_SIZE, 也就是数据量大于网格中所有线程数

const int N = 100000000;
const int M = sizeof(real) * N;
const int BLOCK_SIZE = 128;
const int GRID_SIZE = 10240;//在我们之前的例子中GRID_SIZE直接由(N+BLOCK_SIZE-1)/BLOCK_SIZE计算得到,也就是一个网格算完所有数据

核函数为

__global__ void reduce1(real* d_x, real* d_y, const int N)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int stride = blockDim.x*gridDim.x;//一个网格的所有线程作为步长
    extern __shared__ real s[];
    real temp = 0.0;
    for(int n = bid*blockDim.x+tid; n<N;n+=stride)//在复制之前先进行归约
    {
        temp+=d_x[n];//由于在不同线程块,所以不需要同步
    }
    s[tid]=temp;//根据步长把所有数据归约到gridsize*blocksize长度的数组
    __syncthreads();
    for(int offset = blockDim.x>>1;offset>=32;offset>>=1)
    {
        if(tid<offset)
        {
            s[tid]+=s[tid+offset];
        }
        __syncthreads();
    }
    temp = s[tid];//先把这个数读出来,防止读写竞争
    thread_block_tile<32> g = tiled_partition<32>(this_thread_block());
    //等价于
    // thread_group g32 = tiled_partition(this_thread_block(),32);
    for(int offset = 32; offset>0;offset>>=1)
    {
        temp+=g.shfl_down(temp,offset);
    }//把gridsize*blocksize长度的数组归约到gridsize长度的数组
    if(tid==0)
    {
        d_y[bid]=temp;//0号线程值就是归约和
    }
}

之后调用两次核函数

//把d_x归约到d_y,长度变化为N->GRID_SIZE*BLOCK_SIZE->GRID_SIZE
reduce1<<<GRID_SIZE,BLOCK_SIZE,smem>>>(d_x, d_y,N);
//把d_y归约到d_y,长度变化为GRID_SIZE->1024->1
reduce1<<<1,1024,sizeof(real)*1024>>>(d_y, d_y,GRID_SIZE);

但是当我们循环计算时会频繁动态开辟和释放内存,这占用了大量的时间,一种优化方案是通过静态内存的方式,其在编译时就分配好空间,不会再运行时反复分配,可以如下定义

__device__ real static_y[GRID_SIZE];
real* d_y;
CHECK(cudaGetSymbolAddress((void**)&d_y,static_y));