CUDA性能调优（二）—并行规约及优化

1  并行规约的适用对象

通常用于处理大输入数据集，将一组输入值规约一个值。

数据特点：

（1）对于数据集中的元素没有顺序要求。

（2）可将数据分为若干小集合，每个线程处理一个集合。

操作可以是：求最大值（Max）、求最小值（Min）、求和（Sum）、求乘（Product）。

2  并行规约的实现方法

2.1  简单的求和规约

首先介绍一下最简单的并行规约求和算法。

算法的示意图如下：

      

 对应的程序如下（可直接运行）：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
using namespace std;

const int threadsPerBlock=512 ; 
const int N=2048;
const int blocksPerGrid = (N + threadsPerBlock -1) / threadsPerBlock;

__global__ void ReductionSum(float *d_a, float *d_partial_sum)
{
	//申请共享内存，存在于每个block中 
	__shared__ float partialSum[threadsPerBlock];

	//确定索引
	int i = threadIdx.x + blockIdx.x * blockDim.x;
	int tid = threadIdx.x;

	//传global memory数据到shared memory
	partialSum[tid]=d_a[i];

	//传输同步
	__syncthreads();
	
	//在共享存储器中进行规约
	for(int stride = 1; stride < blockDim.x; stride*=2)
	{
		if(tid%(2*stride)==0) partialSum[tid]+=partialSum[tid+stride];
		__syncthreads();
	}
	//将当前block的计算结果写回输出数组
	if(tid==0)  
		d_partial_sum[blockIdx.x] = partialSum[0];
}

int main()
{
	//申请host端内存及初始化
	float   *h_a,*h_partial_sum;	
    h_a = (float*)malloc( N*sizeof(float) );
    h_partial_sum = (float*)malloc( blocksPerGrid*sizeof(float));
	
	for (int i=0; i < N; ++i)  h_a[i] = i;
	
	//分配显存空间
	int size = sizeof(float);
	float *d_a;
	float *d_partial_sum;
	cudaMalloc((void**)&d_a,N*size);
	cudaMalloc((void**)&d_partial_sum,blocksPerGrid*size);

	//把数据从Host传到Device
	cudaMemcpy(d_a, h_a, size*N, cudaMemcpyHostToDevice);

	//调用内核函数
	ReductionSum<<<blocksPerGrid,threadsPerBlock>>>(d_a,d_partial_sum);

	//将结果传回到主机端
	cudaMemcpy(h_partial_sum, d_partial_sum, size*blocksPerGrid, cudaMemcpyDeviceToHost);

	//将部分和求和
	int sum=0;
    for (int i=0; i < blocksPerGrid; ++i)  sum += h_partial_sum[i];

	cout<<"sum="<<sum<<endl;
	
	//释放显存空间
	cudaFree(d_a);
	cudaFree(d_partial_sum);

	free(h_a);
    free(h_partial_sum);

    return 0;
}

代码分析：

（1）程序中的求模运算效率很低，每次计算至少需要32个时钟周期，如果将求模运算改为等价的乘法指令，则只需要4个指令周期。

（2）上述代码，在一个warp中的线程存在分支，会大大的降低了GPU的执行效率，应尽量避免。

2.2  具有较少线程分支的求和规约

下面这种方案对上面的程序进行了优化，算法示意图如下所示：

 对应的程序如下（可直接运行）：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
using namespace std;

const int threadsPerBlock=512 ; 
const int N=2048;
const int blocksPerGrid = (N + threadsPerBlock -1) / threadsPerBlock;

__global__ void ReductionSum(float *d_a, float *d_partial_sum)
{
	//申请共享内存，存在于每个block中 
	__shared__ float partialSum[threadsPerBlock];

	//确定索引
	int i = threadIdx.x + blockIdx.x * blockDim.x;
	int tid = threadIdx.x;

	//传global memory数据到shared memory
	partialSum[tid]=d_a[i];

	//传输同步
	__syncthreads();
	
	//在共享存储器中进行规约
	for(int stride = blockDim.x/2; stride > 0; stride/=2)
	{
		if(tid<stride) partialSum[tid]+=partialSum[tid+stride];
		__syncthreads();
	}
	//将当前block的计算结果写回输出数组
	if(tid==0)  
		d_partial_sum[blockIdx.x] = partialSum[0];
}

int main()
{
	//申请host端内存及初始化
	float   *h_a,*h_partial_sum;	
    h_a = (float*)malloc( N*sizeof(float) );
    h_partial_sum = (float*)malloc( blocksPerGrid*sizeof(float));
	
	for (int i=0; i < N; ++i)  h_a[i] = 1;
	
	//分配显存空间
	int size = sizeof(float);
	float *d_a;
	float *d_partial_sum;
	cudaMalloc((void**)&d_a,N*size);
	cudaMalloc((void**)&d_partial_sum,blocksPerGrid*size);

	//把数据从Host传到Device
	cudaMemcpy(d_a, h_a, size*N, cudaMemcpyHostToDevice);

	//调用内核函数
	ReductionSum<<<blocksPerGrid,threadsPerBlock>>>(d_a,d_partial_sum);

	//将结果传回到主机端
	cudaMemcpy(h_partial_sum, d_partial_sum, size*blocksPerGrid, cudaMemcpyDeviceToHost);

	//将部分和求和
	int sum=0;
    for (int i=0; i < blocksPerGrid; ++i)  sum += h_partial_sum[i];

	cout<<"sum="<<sum<<endl;
	
	//释放显存空间
	cudaFree(d_a);
	cudaFree(d_partial_sum);

	free(h_a);
    free(h_partial_sum);

    return 0;
}

代码还可继续优化......