FunnyWii
FunnyWii
Published on 2024-06-27 / 68 Visits
0
0

继续学CUDA编程(三) —— CUDA执行模型(下)

并行归约问题

先解释一下什么是归约(Reduction),归约是将某个计算问题变换为另一个问题的过程。在CUDA运算中,在向量中执行满足交换律和结合律的运算,被称为归约问题。每次迭代计算方式都是相同的(归),从一组多个数据最后得到一个数(约)^[1]^。比如当给定N个数值,求其SUM/MAX/MIN/MEAN 值的操作,也可被称为归约。

如果在CPU上,我们只需要一个循环遍历就可以实现以上操作,但使用GPU呢?

  1. 将输入划分到更小的数据块中
  2. 用一个Thread计算一个数据块的部分和
  3. 对每个数据块的部分和,再求和,得出最终结果

上述方法通常使用迭代成对实现,具体过程见下图:

并行规约.png

可以看出,用了三轮迭代就计算出了数组总和。

并行归约中的分化

上图是使用相邻配对方法实现的,这里有两个全局内存数组,一个大数组存放整个数组并进行归约;另一个小数组用于存放每个Block的部分和。由于归约是就地完成的,意味着在每一步全局内存中都有部分值被替换。使用 _syncthreads()可以保证Block中的任一线程在进入下次迭代之前,当前迭代中的所有部分和都被保存在全局内存中。

//这个g_idata存放的是整个数组,进行存放整个数组;g_odata 全局内存中小数组,存放每个线程块的部分和
__global__ void reduceNeighbored(int* g_idata, int* g_odata, unsigned int n) {
	//首先获取thread ID
	unsigned int tid = threadIdx.x;
	//获取idx全局编号
	unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;

	//将全局的数据指针转化为这个块的局部数据指针 idata指向某个块的开始
	int* idata = g_idata + blockIdx.x * blockDim.x;

	//数组越界检查
	if (idx >= n) {
		return;
	}

	//就地对全局数据进行规约
	for (int stride = 1; stride < blockDim.x; stride *= 2) {
		if ((tid % (2 * stride)) == 0) {
			idata[tid] += idata[tid + stride];
		}
		//它能保证线程块之间的同步问题,就是线程块中的任一线程在进入下一步迭代之前,在当前迭代里每个线程的所有部分和都被保存在全局内存
		__syncthreads();
	}
    // Block之间不能同步,因此要将Block产生的部分和返回Host
	// 将块的局部结果写到全局内存  blockIdx的范围[0,grid.x)
	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];
}

其中 (tid % (2 * stride)) == 0这个判断可能不太好理解,解释一下:假设和上图一样,8个数字相加。tid编号为0-7,最开始stride=1,乘2后,只有0、2、4、6四个tid满足条件,因此这4个位置的内存,操作 idata[tid] = idata[tid] + idata[tid+stride];随后,stride \times 2=2,只有0、2满足条件,同样执行加法操作,最后,只有0满足条件。

让我们细化到第一轮:
tid 为 0 的线程:idata[0] += idata[1](有效,因为 0 % 2 == 0)
tid 为 1 的线程:不执行求和操作(因为 1 % 2 != 0)
tid 为 2 的线程:idata[2] += idata[3](有效,因为 2 % 2 == 0)
tid 为 3 的线程:不执行求和操作(因为 3 % 2 != 0)
tid 为 4 的线程:idata[4] += idata[5](有效,因为 4 % 2 == 0)
tid 为 5 的线程:不执行求和操作(因为 5 % 2 != 0)
tid 为 6 的线程:idata[6] += idata[7](有效,因为 6 % 2 == 0)
tid 为 7 的线程:不执行求和操作(因为 7 % 2 != 0)

改善分化

由于 if ((tid % (2 * stride)) == 0) 只对偶数ID线程为 true,即只有偶数ID线程会被调度,其他线程闲置,因此会导致线程束分化。可以将访问索引修改为 int index = 2 ∗ stride ∗ tid 。注意这里仍是相邻配对,而非交错配对。

当跨度乘2后,index可以后移到有数据的内存,而非使用tid对应内存地址,减少了线程空闲。首先我们保证在一个块中前几个执行的线程束是在接近满跑的,而后半部分线程束基本是不需要执行的,当一个线程束内存在分支,而分支都不需要执行的时候,硬件会停止他们调用别人,这样就节省了资源,所以高效体现在这,如果还是所有分支不满足的也要执行,即便整个线程束都不需要执行的时候,那这种方案就无效了^[1]^。

让我们细化到第一轮:
tid 为 0 的线程:idata[0] += idata[1](有效,因为 0 < 8)
tid 为 1 的线程:idata[2] += idata[3](有效,因为 2 < 8)
tid 为 2 的线程:idata[4] += idata[5](有效,因为 4 < 8)
tid 为 3 的线程:idata[6] += idata[7](有效,因为 6 < 8)
tid 为 4 的线程:尝试 idata[8] += idata[9](无效,因为 8 >= 8)
tid 为 5 的线程:尝试 idata[10] += idata[11](无效,因为 10 >= 8)
tid 为 6 的线程:尝试 idata[12] += idata[13](无效,因为 12 >= 8)
tid 为 7 的线程:尝试 idata[14] += idata[15](无效,因为 14 >= 8)

通过对第一轮操作的细化,应该能够理解,为什么 if ((tid % (2 * stride)) == 0) 可以让Block的前半部分执行求和操作。

交错配对归约

交错配对归约.png

交错配对的方法,初始stride是Block大小的一半,然后每次迭代减少一半。

核心部分代码如下:

for(int stride = blockDim.x /2; stride >0; stride >>1){
	if(tid < stride){
    	idata[tid] += idata[tid] + idata[tid+stride]
    }
}

其中 stride >>1,它将 stride 的二进制表示向右移动 1 位。移出的位被丢弃,假设初始值为8,一次循环后则变为4。

展开循环

封闭的循环,可以将其迭代次数减少甚至完全消除。循环体的复制次数为循环展开因子,迭代次数变为原迭代次数/循环展开因子

一个比较简单的例子:

for (int i = 0; i < 100; i++){
	a[i] = b[i] + c[i];
}
for (int i = 0; i < 50; i +=2){
	a[i] = b[i] + c[i];
	a[i+1] = b[i+1] + c[i+1];
}

展开的归约

在之前的归约函数中,每个Block只能处理一部分数据,但是如果能用一个Block处理更多数据,即在数据求和前,使用每个Thread进行一次加法,再从别的Block取数据,相当于先做一次向量加法,再进行归约。

__global__ void reduceUnroll2(int * g_idata,int * g_odata,unsigned int n){
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x*2+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x*2;
	if(idx+blockDim.x<n){
		g_idata[idx]+=g_idata[idx+blockDim.x];
	}
	__syncthreads();
	//in-place reduction in global memory
	for (int stride = blockDim.x/2; stride>0 ; stride >>=1){
		if (tid <stride){
			idata[tid] += idata[tid + stride];
		}
		//synchronize within block
		__syncthreads();
	}
	//write result for this block to global mem
	if (tid == 0)  g_odata[blockIdx.x] = idata[0];
}

g_idata[idx]+=g_idata[idx+blockDim.x]; 这行代码的作用是让每个Thread都与相邻数据块同位置的元素相加。

代码中 unsigned int idx = blockDim.x*blockIdx.x*2+threadIdx.x; 以及 int *idata = g_idata + blockIdx.x*blockDim.x*2; 在确定Block对应数据位置时有个乘2的偏移量。

展开的归约.png

根据上图,因为只需要一半的Block来处理数据,因此全局数据索引也被调整

展开线程的归约

当只剩下32个或更少线程时,由于Warp执行是SIMT,每条指令后有隐式的Warp内同步过程。归约循环的最后6个迭代可以用下面语句:

if(tid<32){
		volatile int *vmem = idata;
		vmem[tid]+=vmem[tid+32];
		vmem[tid]+=vmem[tid+16];
		vmem[tid]+=vmem[tid+8];
		vmem[tid]+=vmem[tid+4];
		vmem[tid]+=vmem[tid+2];
		vmem[tid]+=vmem[tid+1];
	}

这个Warp展开避免了执行循环控制和线程同步。注意:一个Warp的执行进度是完全相同的,当执行 vem[tid+32] 的时候,这32个线程都在执行这步。本Warp内的所有线程都不会执行到 vem[tid+16]

所以要添加代码中的 volatile ,这个修饰符用来修饰变量 vmem,让编译器每次赋值时必须将 vmem[tid]的值存回全局内存中。控制变量结果写回到内存,而不是存在共享内存,或者缓存中,因为下一步的计算马上要用到它,如果写入缓存,可能造成下一步的读取会读到错误的数据。

当只剩下最后的64个数合并成1个数时,首先将前32个数,按照stride=32来进行并行加法,得到64个数字的两两之和,存在前32个数字中。

然后这32个数加上步长为16的变量,理论上,这样能得到16个数,这16个数的和就是最后这个块的归约结果^[2]^。

完全展开的归约

如果已知迭代次数,那么可以把循环完全展开。

__global__ void reduceUnrollWarp8(int * g_idata,int * g_odata,unsigned int n){
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x*8+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x*8;
	//unrolling 8;
	if(idx+7 * blockDim.x<n){
		int a1=g_idata[idx];
		int a2=g_idata[idx+blockDim.x];
		int a3=g_idata[idx+2*blockDim.x];
		int a4=g_idata[idx+3*blockDim.x];
		int a5=g_idata[idx+4*blockDim.x];
		int a6=g_idata[idx+5*blockDim.x];
		int a7=g_idata[idx+6*blockDim.x];
		int a8=g_idata[idx+7*blockDim.x];
		g_idata[idx]=a1+a2+a3+a4+a5+a6+a7+a8;
	}
	__syncthreads();
	//in-place reduction in global memory
	for (int stride = blockDim.x/2; stride>32; stride >>=1){
		if (tid <stride){
			idata[tid] += idata[tid + stride];
		}
		//synchronize within block
		__syncthreads();
	}
	//write result for this block to global mem
	if(tid<32){
		volatile int *vmem = idata;
		vmem[tid]+=vmem[tid+32];
		vmem[tid]+=vmem[tid+16];
		vmem[tid]+=vmem[tid+8];
		vmem[tid]+=vmem[tid+4];
		vmem[tid]+=vmem[tid+2];
		vmem[tid]+=vmem[tid+1];
	}

	if (tid == 0)  g_odata[blockIdx.x] = idata[0];
}

模板函数的归约

既然展开了最后64个数据,那也可以展开更多数据。在 __device__函数中,CUDA支持模板参数,我们可以指定展开的数量作为参数。

template <unsigned int iBlockSize>

if(iBlockSize>=1024 && tid <512)
    idata[tid]+=idata[tid+512];
__syncthreads();
if(iBlockSize>=512 && tid <256)
    idata[tid]+=idata[tid+256];
__syncthreads();
if(iBlockSize>=256 && tid <128)
    idata[tid]+=idata[tid+128];
__syncthreads();
if(iBlockSize>=128 && tid <64)
    idata[tid]+=idata[tid+64];
__syncthreads();

对于这部分代码,编译器会去检查 blockDim.x 是否固定。如果固定,则会移除不可能实现的部分。而且当 blockDim.x 确定后就不能更改了。

参考文章

[1] 【CUDA 基础】3.4 避免分支分化
[2] CUDA编程权威指针(13)-改善并行规约的分化


Comment