高性能计算-GPU并行规约(27)
1. 目标：对数组进行求和，并做优化对比
2. baseline 代码
相邻求和: 根据blockSize对数据分块,并将数据放在共享内存,以线程块为单位,块内线程数量=数据个数,相邻配对,用其中第一个元素索引为ID的线程进行计算,计算结果放在第一个元素位置,循环进行下一轮计算,最后块求和计算结果赋值到全局内存以blockIdx.x为id的数组中,拷贝到主机对数组求和得到所有数字之和。
#include <stdio.h>
#include "common.h"
#define BLOCKSIZE 512
__global__ void kernel_reduce(float* in,long N,float* out,int gridSize)
{
//每个线程块加载与线程数量相同的数据
__shared__ float arrShared[BLOCKSIZE];
int tid = threadIdx.x;
//共享内存数据初始化
arrShared[tid] = in[blockIdx.x*blockDim.x + tid];
__syncthreads();
//线程块内计算
for(long s=1;s<=BLOCKSIZE/2;s*=2)
{
if(0 == tid % (2*s))
arrShared[tid] += arrShared[tid+s];
__syncthreads();
}
if(tid ==0)
out[blockIdx.x] = arrShared[0];
}
int main(int argc, char ** argv)
{
long N =1<<10;
if(argc > 1)
N = 1<<atoi(argv[1]);
N *= 32;
//GPU计算参数
float gpuTime = 0;
int gridSize = (N-1)/BLOCKSIZE+1; //数组初始化
printf("N %ld ,GridSize %d\n",N,gridSize);
//为了保证对比准确性,最后求和不计入耗时对比
float cpuTime = 0;
float* arrHost = (float*)calloc(N,sizeof(float));
float* arrResult = (float*)calloc(gridSize,sizeof(float));
float resultHost = 0;
initialData(arrHost,N);
double start = cpuSecond();
for(int i=0;i<N/BLOCKSIZE;i++)
{
float temp =0;
for(int j=0;j<BLOCKSIZE;j++)
temp += arrHost[i*BLOCKSIZE+j];
arrResult[i] = temp;
}
for(int i=0;i<N%BLOCKSIZE;i++)
arrResult[gridSize-1] += arrHost[N/BLOCKSIZE*BLOCKSIZE + i];
double end = cpuSecond();
cpuTime = (end - start)*1000;
//gpu
float *arrD = NULL;
float *resultD = NULL;
float *resultFromD = NULL;
float resultGpu = 0;
CudaSafeCall(cudaMalloc((void**)&arrD,N*sizeof(float)));
CudaSafeCall(cudaMalloc((void**)&resultD,gridSize*sizeof(float)));
resultFromD = (float*)calloc(gridSize,sizeof(float));
cudaEvent_t startD;
cudaEvent_t endD;
CudaSafeCall(cudaEventCreate(&startD));
CudaSafeCall(cudaEventCreate(&endD));
CudaSafeCall(cudaEventRecord(startD));
CudaSafeCall(cudaMemcpy(arrD,arrHost,N*sizeof(float),cudaMemcpyHostToDevice));
kernel_reduce<<<gridSize,BLOCKSIZE,sizeof(float)*BLOCKSIZE>>>(arrD,N,resultD,gridSize);
CudaCheckError();
CudaSafeCall(cudaMemcpy(resultFromD,resultD,gridSize*sizeof(float),cudaMemcpyDeviceToHost));
CudaSafeCall(cudaEventRecord(endD));
cudaEventSynchronize(endD);
CudaSafeCall(cudaEventElapsedTime(&gpuTime,startD,endD));
CudaSafeCall(cudaEventDestroy(startD));
CudaSafeCall(cudaEventDestroy(endD));
//汇总求和
for(int i=0;i<gridSize;i++)
{
resultHost += arrResult[i];
resultGpu += resultFromD[i];
}
printf("数据量 %ld ;串行结算结果为%.3f,耗时 %.3f ms;GPU计算结果为%.3f,耗时 %.3f ms;加速比为%.3f\n",N,resultHost,cpuTime,resultGpu,gpuTime,cpuTime/gpuTime);
CudaSafeCall(cudaFree(arrD));
CudaSafeCall(cudaFree(resultD));
free(arrHost);
free(resultFromD);
return 0;
}
3. 优化代码及思路
优化一:从全局内存加载时计算,一个线程块计算8个线程块的数据
优化二:从全局内存加载数据时使用合并访问,一个线程处理的数据索引步长大小为 blocksize 大小,
这样处理的话在一个线程束中访存可以对 cacheline 连续命中
优化三:线程束内做循环展开
优化四:用 shfl API 线程束内访问其他线程寄存器
/*
数组求和
优化一:从全局内存加载时计算,一个线程块计算8个线程块的数据
优化二:从全局内存加载数据时使用合并访问,一个线程处理的数据索引步长大小为 blocksize 大小,
这样处理的话在一个线程束中访存可以对 cacheline 连续命中
优化三:线程束内做循环展开
优化四:用 shfl API 线程束内访问其他线程寄存器
注意：如果是 float 每个block cpu和gpu计算数据顺序不同，大数+小数,小数可能会被吃掉,精度不能太高
*/
#include <stdio.h>
#include "common.h"
#define BLOCK_SIZE 256 // 最大1024
#define NPerThread 8 // 每个线程初始化数据个数
#define WARP_SIZE 32 // 线程束大小
// #define warpNum 1024/WARP_SIZE
// 优化三:线程束内的规约,并做循环展开
// 优化四:使用 shfl API 可以在warp内跨线程访问寄存器数据
template <int warp_len> // warp 中有效数据个数
__device__ float warp_reduce(float sum)
// __device__ float warpReduceSum(float sum)
{
// 防止线程块大小设置太小,做判断
#if 0
if (warp_len >= 32)
// 向前面的线程传递寄存器数值
sum += __shfl_down_sync(0xffffffff, sum, 16);
if (warp_len >= 16)
sum += __shfl_down_sync(0xffffffff, sum, 8);
if (warp_len >= 8)
sum += __shfl_down_sync(0xffffffff, sum, 4);
if (warp_len >= 4)
sum += __shfl_down_sync(0xffffffff, sum, 2);
if (warp_len >= 2)
sum += __shfl_down_sync(0xffffffff, sum, 1);
#else
for (int i = warp_len >> 1; i >= 1; i >>= 1)
sum += __shfl_down_sync(0xFFFFFFFF, sum, i);
#endif
return sum;
}
#if 1
// 线程块大小 每个线程要处理的个数
template <int blockSize, int NUM_PER_THREAD>
__global__ void block_reduce(float *g_in, float *g_out)
{
float sum = 0; // 保存当前线程加和的数值
int tid = threadIdx.x;
// 两种访存方式都没有 bankconflict ，共享内存访存事务数量一致
// 写法二全局内存非合并访存 wavefront 比 写法一增加，但是L1缓存命中高，意义不大，因为共享内存访存事务数量一致
#if 1 // 写法一
// 优化二:全局内存合并访存，无 bankconflict
// 在 NUM_PER_THREAD 个线程块相同位置分别取数据给线程求和
// 而不是对一个块中连续内存访问求和
int temp1 = blockIdx.x * blockSize * NUM_PER_THREAD + tid;
#pragma unrool
for (int i = 0; i < NUM_PER_THREAD; i++)
// 优化一:加载时计算
sum += g_in[temp1 + i * blockSize];
#else // 写法二
// 优化: 全局内存非合并访存，无 bankconflict
// 从全局内存取数据时,每个线程取多个数据
int tempId1 = blockIdx.x * blockDim.x * NUM_PER_THREAD;
int tempId2 = tid * NUM_PER_THREAD;
// 共享内存数据初始化,循环展开
#pragma unrool
for (int i = 0; i < NUM_PER_THREAD; i++)
sum += g_in[tempId1 + tempId2 + i];
#endif
// 共享内存默认初始化为0
// __shared__ float arrWarp[WARP_SIZE];
__shared__ float arrWarp[blockSize / WARP_SIZE];
// 做warp内规约
// 使用 shfl 接口,传入要从其他线程复制的寄存器变量 sum;每个线程都要参与交换数据
sum = warp_reduce<WARP_SIZE>(sum);
// wrap内 0号线程将结果保存到共享内存
if (tid % WARP_SIZE == 0)
arrWarp[tid / WARP_SIZE] = sum;
__syncthreads();
// 将warp规约结果放入第一个warp内
sum = (tid < blockSize / WARP_SIZE) ? arrWarp[tid] : 0;
// 对第一个warp规约（共享内存规约,blockDim.x最大1024，1024/32=32）
if ((tid / WARP_SIZE) == 0)
sum = warp_reduce<blockSize / WARP_SIZE>(sum);
// 保存整个线程块规约结果
if (tid == 0)
g_out[blockIdx.x] = sum;
}
#else
template <int blockSize, int NUM_PER_THREAD>
__global__ void block_reduce(float *g_in, float *g_out)
{
int tid = threadIdx.x;
float sum = 0.0;
// 共享内存初始化
// 一个线程获取NUM_PER_THREAD个数据求和
int offset = blockIdx.x * blockDim.x * NUM_PER_THREAD;
for (int i = 0; i < NUM_PER_THREAD; i++)
sum += g_in[offset + i * blockDim.x + tid];
// 计算
// warp内规约
__shared__ float arrWarp[blockSize / WARP_SIZE];
sum = warp_reduce<32>(sum);
if (tid % warpSize == 0)
arrWarp[tid / warpSize] = sum;
__syncthreads();
// 共享内存规约,blockDim.x最大1024，1024/32=32
if (tid < blockDim.x / warpSize)
sum = arrWarp[tid];
else
sum = 0;
// 共享内存最大32个元素，只需用第一个warp做规约
if (tid / WARP_SIZE == 0)
sum = warp_reduce<blockSize / WARP_SIZE>(sum);
// 写回数据
if (tid == 0)
g_out[blockIdx.x] = sum;
}
#endif
int main(int argc, char **argv)
{
// printf("0\n");
long int N = 1 << 15;
if (argc > 1)
N = 1 << atoi(argv[1]);
// GPU参数计算
float gpuTime = 0;
// printf("1\n");
int gridSize = (N - 1) / (BLOCK_SIZE * NPerThread) + 1;
// int gridSize = (N-1)/BLOCK_SIZE + 1;
printf("N %ld ,GridSize %d\n", N, gridSize);
// cpu
float cpuTime = 0;
float *arrHost = (float *)calloc(N, sizeof(float));
float *arrResult = (float *)calloc(gridSize, sizeof(float));
float resultHost = 0;
initialData(arrHost, N);
double start = cpuSecond();
int number = BLOCK_SIZE * NPerThread;
for (int i = 0; i < N / number; i++)
{
float temp = 0;
for (int j = 0; j < number; j++)
temp += arrHost[i * number + j];
arrResult[i] = temp;
}
for (int i = 0; i < N % (BLOCK_SIZE * NPerThread); i++)
arrResult[gridSize - 1] += arrHost[N / number * number + i];
double end = cpuSecond();
cpuTime = (end - start) * 1000;
// gpu
float *arrD = NULL;
float *resultD = NULL;
float *resultFromD = NULL;
float resultGpu = 0;
CudaSafeCall(cudaMalloc((void **)&arrD, N * sizeof(float)));
CudaSafeCall(cudaMalloc((void **)&resultD, gridSize * sizeof(float)));
resultFromD = (float *)calloc(gridSize, sizeof(float));
cudaEvent_t startD;
cudaEvent_t endD;
CudaSafeCall(cudaEventCreate(&startD));
CudaSafeCall(cudaEventCreate(&endD));
CudaSafeCall(cudaMemcpy(arrD, arrHost, N * sizeof(float), cudaMemcpyHostToDevice));
CudaSafeCall(cudaEventRecord(startD, 0));
// for (int i = 0; i < 6; i++)
block_reduce<BLOCK_SIZE, NPerThread><<<gridSize, BLOCK_SIZE, BLOCK_SIZE / WARP_SIZE * sizeof(float)>>>(arrD, resultD);
CudaCheckError();
CudaSafeCall(cudaEventRecord(endD, 0));
CudaSafeCall(cudaMemcpy(resultFromD, resultD, gridSize * sizeof(float), cudaMemcpyDeviceToHost));
cudaEventSynchronize(endD);
CudaSafeCall(cudaEventElapsedTime(&gpuTime, startD, endD));
CudaSafeCall(cudaEventDestroy(startD));
CudaSafeCall(cudaEventDestroy(endD));
// 汇总求和
// 如果是 float 每个block cpu和gpu计算数据顺序不同，大数+小数,小数可能会被吃掉,精度不能太高
for (int i = 0; i < gridSize; i++)
{
resultHost += arrResult[i];
resultGpu += resultFromD[i];
}
printf("数据量 %ld ;串行结算结果为%.3f,耗时 %.3f ms;GPU计算结果为%.3f,耗时 %.3f ms;加速比为%.3f\n", N, resultHost, cpuTime, resultGpu, gpuTime, cpuTime / gpuTime);
CudaSafeCall(cudaFree(arrD));
CudaSafeCall(cudaFree(resultD));
free(arrHost);
free(arrResult);
free(resultFromD);
return 0;
}
5. 测试结果
| 串行(ms) | gpu(ms) | 加速比 | |
|---|---|---|---|
| baseline | 95.2 | 5.3 | 17.8 |
| 优化 | 89.1 | 1.029 | 87 |
| thrust | 1.158 |
6. 结果分析
（1）baseline版本相邻求和，时间复杂度从 O(n)，降为log(n)。
（2）经测试优化中的优化二合并访问内存能加速3倍，还有加载时计算也对性能提升较大。
（3）经过优化后的代码达到了与 thrust reduce差不多的效率。
7. 注意
（1）如果使用float 数据类型规约，CPU和GPU对一个block处理的数据顺序不同，大数加小数会发生大数吃小数的情况，计算精度不一样，误差精度不应该设置太高。
浙公网安备 33010602011771号