CUDA并行规约(交错配对-展开规约）

之前的规约计算，每个线程块负责一个对应的数据块。现在将每个线程块负责两个数据块的规约，从而达到消除指令消耗，增加更多的独立指令的调度，来提高性能。下面是展开因子为2的示意图，将相邻的数据块加到当前线程块对应的数据块，然后进行规约求和。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include "math.h"
#include "stdlib.h"

//错误检查的宏定义
#define CHECK(call)									\
{													\
	const cudaError_t status=call;					\
if (status!=cudaSuccess)							\
	{												\
	printf("文件:%s,函数:%s,行号:%d\n",__FILE__,		\
						__FUNCTION__,__LINE__);		\
	printf("%s", cudaGetErrorString(status));		\
	exit(1);										\
	}												\
}													\

//核函数
__global__ void Kernel(int *d_data, int *d_local_sum, int N)
{
	int tid = threadIdx.x;
	int index = 2*blockIdx.x*blockDim.x + threadIdx.x;
	int *data = d_data +2*blockIdx.x*blockDim.x;


	if (index + blockDim.x < N) d_data[index] += d_data[index+blockDim.x];

	__syncthreads();

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

		__syncthreads();
	}


	if (tid == 0)
	{
		d_local_sum[blockIdx.x] = data[0];
	}
}

//主函数
int main()
{

	//基本参数设置
	cudaSetDevice(0);
	const int N = 1 << 29;
	int local_length = 1024;
	long long total_sum = 0;

	dim3 grid(((N + local_length - 1) / local_length), 1);
	dim3 block(local_length, 1);

	int *h_data = nullptr;
	int *h_local_sum = nullptr;
	int *d_data = nullptr;
	int *d_local_sum = nullptr;


	//Host&Deivce内存申请及数组初始化
	h_data = (int*)malloc(N * sizeof(int));
	h_local_sum = (int*)malloc(int(grid.x/2) * sizeof(int));


	CHECK(cudaMalloc((void**)&d_data, N * sizeof(int)));

	CHECK(cudaMalloc((void**)&d_local_sum, int(grid.x / 2) * sizeof(int)));

	for (int i = 0; i < N; i++)
		h_data[i] = int(10 * sin(0.02*3.14*i));//限制数组元素值，防止最终求和值超过long long的范围

	//数据拷贝至Device
	CHECK(cudaMemcpy(d_data, h_data, N * sizeof(int), cudaMemcpyHostToDevice));

	//for (int i = 0; i < 200; i++)
		//执行核函数
		Kernel << <grid.x/2, block >> > (d_data, d_local_sum, N);

	//数据拷贝至Host
	CHECK(cudaMemcpy(h_local_sum, d_local_sum, int(grid.x/2) * sizeof(int),
		cudaMemcpyDeviceToHost));

	//同步&重置设备
	CHECK(cudaDeviceSynchronize());
	CHECK(cudaDeviceReset());


	for (int i = 0; i < int(grid.x / 2); i++)
	{
		total_sum += h_local_sum[i];
	}

	printf("%I64d \n", total_sum);

	//getchar();
	return 0;

}

同理，可以将展开因子调整为4和8，下面给出展开因子为4的示意图，8的线条太多了。

展开因子为4的代码：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include "math.h"
#include "stdlib.h"

//错误检查的宏定义
#define CHECK(call)									\
{													\
	const cudaError_t status=call;					\
if (status!=cudaSuccess)							\
	{												\
	printf("文件:%s,函数:%s,行号:%d\n",__FILE__,		\
						__FUNCTION__,__LINE__);		\
	printf("%s", cudaGetErrorString(status));		\
	exit(1);										\
	}												\
}													\

//核函数
__global__ void Kernel(int *d_data, int *d_local_sum, int N)
{
	int tid = threadIdx.x;
	int index = 4 * blockIdx.x*blockDim.x + threadIdx.x;
	int *data = d_data + 4 * blockIdx.x*blockDim.x;


	if (index + 3 * blockDim.x < N)
	{
		int a = d_data[index];
		int a1 = d_data[index + blockDim.x];
		int a2 = d_data[index + 2 * blockDim.x];
		int a3 = d_data[index + 3 * blockDim.x];
		d_data[index] = (a+ a1+ a2 +a3 );
	}

	__syncthreads();

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

		__syncthreads();
	}


	if (tid == 0)
	{
		d_local_sum[blockIdx.x] = data[0];
	}
}

//主函数
int main()
{

	//基本参数设置
	cudaSetDevice(0);
	const int N = 1 << 29;
	int local_length = 1024;
	long long total_sum = 0;

	dim3 grid(((N + local_length - 1) / local_length), 1);
	dim3 block(local_length, 1);

	int *h_data = nullptr;
	int *h_local_sum = nullptr;
	int *d_data = nullptr;
	int *d_local_sum = nullptr;


	//Host&Deivce内存申请及数组初始化
	h_data = (int*)malloc(N * sizeof(int));
	h_local_sum = (int*)malloc(int(grid.x / 4) * sizeof(int));


	CHECK(cudaMalloc((void**)&d_data, N * sizeof(int)));

	CHECK(cudaMalloc((void**)&d_local_sum, int(grid.x / 4) * sizeof(int)));

	int s = 0;
	for (int i = 0; i < N; i++)
	{
		h_data[i] = int(10 * sin(0.02*3.14*i));//限制数组元素值，防止最终求和值超过long long的范围
		s += h_data[i];
	}

											   //数据拷贝至Device
	CHECK(cudaMemcpy(d_data, h_data, N * sizeof(int), cudaMemcpyHostToDevice));

	//for (int i = 0; i < 200; i++)
	//执行核函数
	Kernel << <grid.x / 4, block >> > (d_data, d_local_sum, N);

	//数据拷贝至Host
	CHECK(cudaMemcpy(h_local_sum, d_local_sum, int(grid.x / 4) * sizeof(int),
		cudaMemcpyDeviceToHost));

	//同步&重置设备
	CHECK(cudaDeviceSynchronize());
	CHECK(cudaDeviceReset());


	for (int i = 0; i < int(grid.x/4); i++)
	{
		total_sum += h_local_sum[i];
	}

	printf("%I64d \n", total_sum);

	//getchar();
	return 0;

}

展开因子为8的代码：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include "math.h"
#include "stdlib.h"

//错误检查的宏定义
#define CHECK(call)									\
{													\
	const cudaError_t status=call;					\
if (status!=cudaSuccess)							\
	{												\
	printf("文件:%s,函数:%s,行号:%d\n",__FILE__,		\
						__FUNCTION__,__LINE__);		\
	printf("%s", cudaGetErrorString(status));		\
	exit(1);										\
	}												\
}													\

//核函数
__global__ void Kernel(int *d_data, int *d_local_sum, int N)
{
	int tid = threadIdx.x;
	int index = 8 * blockIdx.x*blockDim.x + threadIdx.x;
	int *data = d_data + 8 * blockIdx.x*blockDim.x;


	if (index +	7 * blockDim.x < N)
	{
		int a = d_data[index];
		int a1 = d_data[index + blockDim.x];
		int a2 = d_data[index + 2 * blockDim.x];
		int a3 = d_data[index + 3 * blockDim.x];
		int a4 = d_data[index + 4 * blockDim.x];
		int a5 = d_data[index + 5 * blockDim.x];
		int a6 = d_data[index + 6 * blockDim.x];
		int a7 = d_data[index + 7 * blockDim.x];
		d_data[index] = (a + a1 + a2 + a3+ a4+ a5+ a6 + a7);
	}

	__syncthreads();

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

		__syncthreads();
	}


	if (tid == 0)
	{
		d_local_sum[blockIdx.x] = data[0];
	}
}

//主函数
int main()
{

	//基本参数设置
	cudaSetDevice(0);
	const int N = 1 << 29;
	int local_length = 1024;
	long long total_sum = 0;

	dim3 grid(((N + local_length - 1) / local_length), 1);
	dim3 block(local_length, 1);

	int *h_data = nullptr;
	int *h_local_sum = nullptr;
	int *d_data = nullptr;
	int *d_local_sum = nullptr;


	//Host&Deivce内存申请及数组初始化
	h_data = (int*)malloc(N * sizeof(int));
	h_local_sum = (int*)malloc(int(grid.x / 8) * sizeof(int));


	CHECK(cudaMalloc((void**)&d_data, N * sizeof(int)));

	CHECK(cudaMalloc((void**)&d_local_sum, int(grid.x / 8) * sizeof(int)));

	int s = 0;
	for (int i = 0; i < N; i++)
	{
		h_data[i] = int(10 * sin(0.02*3.14*i));//限制数组元素值，防止最终求和值超过long long的范围
		s += h_data[i];
	}

	//数据拷贝至Device
	CHECK(cudaMemcpy(d_data, h_data, N * sizeof(int), cudaMemcpyHostToDevice));

	//for (int i = 0; i < 200; i++)
	//执行核函数
	Kernel << <grid.x / 8, block >> > (d_data, d_local_sum, N);

	//数据拷贝至Host
	CHECK(cudaMemcpy(h_local_sum, d_local_sum, int(grid.x / 8) * sizeof(int),
		cudaMemcpyDeviceToHost));

	//同步&重置设备
	CHECK(cudaDeviceSynchronize());
	CHECK(cudaDeviceReset());


	for (int i = 0; i < int(grid.x / 8); i++)
	{
		total_sum += h_local_sum[i];
	}

	printf("%I64d \n", total_sum);

	//getchar();
	return 0;

}