高性能计算-CUDA一维信号均值滤波及内存优化(22)
1. 目标：使用CPU和GPU对一千万数量级的一维信号进行均值滤波，并且根据GPU存储模型对数据存储进行优化，最终对比计算结果并计算加速比。
2. 代码
/*
cuda实现对一维信号卷积平滑滤波处理，并于串行计算对比结果和加速比,卷积核大小为5
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifdef WIN32
#include <sys/utime.h>
#else
#include <sys/time.h>
#endif
#include "cuda_runtime_api.h"
#define MAXLEN 256
#define QW 10000000 //数据量
#define EP 1e-15 //精度
#define FSIZE 5 //一维卷积核大小
#define BlockSize 512
double* in = (double*)calloc(QW,sizeof(double));
double* out = (double*)calloc(QW, sizeof(double)); //cpu串行计算结果
double* outFD = (double*)calloc(QW, sizeof(double)); //gpu计算结果
double* filter = (double*)calloc(FSIZE, sizeof(double));
__constant__ double fliterConst[FSIZE];
void readFile(double*in,int len)
{
FILE* file = fopen("./noise1qw.txt", "r");
if (!file)
{
perror("打开noise失败");
return;
}
//fgets读取一行包含换行符和字符串结束字符\0，这里使用fscanf读取遇到空格和换行符结束
int i = 0;
while (fscanf(file,"%lf",in+i++) >0)
{
//printf("%1.16f\n", in[i - 1]);
//if (i == 2)
//break;
}
fclose(file);
}
void saveFile(double* out, int len)
{
FILE* file = fopen("./result.txt", "w");
if (file == NULL)
{
perror("open result.txt failed.\n");
return;
}
for (int i = 0;i < len;i++)
fprintf(file, "%2.16f\n", out[i]);
fclose(file);
}
//ms
long getTime()
{
struct timeval cur;
gettimeofday(&cur, NULL);
// printf("sec %ld usec %ld,toal ms %ld\n",cur.tv_sec,cur.tv_usec,cur.tv_sec*1e3 + cur.tv_usec / 1e3);
return cur.tv_sec*1e3 + cur.tv_usec / 1e3;
}
//串行卷积
void serialConv(double* in,double* out,size_t len,double* filter,size_t filterSize)
{
for (int i = 0;i < len;i++)
{
int start = i - FSIZE / 2;
double result = 0.0;
for (int j = 0;j < FSIZE;j++)
{
int offPos = start + j;
if (offPos >= 0 && offPos < len)
result += filter[j] * in[offPos];
}
out[i] = result;
}
}
//gpu并行卷积
__global__ void kernalConv(double* in,double* out, size_t len, double* filter, size_t filterSize)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
int start = id - FSIZE / 2;
double result = 0.0;
for (int i = 0;i < FSIZE;i++)
{
int offPos = start + i;
if (offPos >= 0 && offPos < len)
result += filter[i] * in[offPos];
}
out[id] = result;
}
//优化1：常用卷积核常量
__global__ void kernalConv_const(double* in,double* out, size_t len)
{
int id = blockIdx.x * blockDim.x + threadIdx.x;
int start = id - FSIZE / 2;
double result = 0.0;
for (int i = 0;i < FSIZE;i++)
{
int offPos = start + i;
if (offPos >= 0 && offPos < len)
result += fliterConst[i] * in[offPos];
}
out[id] = result;
}
//优化1：使用卷积核放常量
//优化2：将一个线程块中可能访问的输入数据放入共享内存，防止频繁访问全局内存
//分析：每个共享内存变量中存放本block需要计算的数据
__global__ void kernelConv_const_share(double* in,double* out, size_t len)
{
//每个block中都有此shareM数据空间，需要每个块分别初始化自己的计算数据
//共享内存初始化思路：在块内，根据threadIdx判断当前线程在块内是否前/后raduis个线程，需要额外初始化padding数据；
// (1)如果是索引为0 的块，左侧的padding应为0，否则应为前一个数据块的后radius个数据
// (2)如果是索引为最后一个块，右侧的padding应为0，否则应为后一个数据块的前radius个数据
//并且每个线程根据全局唯一一维id作为输入数据索引，直接数据赋值给共享内存[threadIdx+radius]索引的空间
__shared__ double shareM[BlockSize+FSIZE-1];
int id = blockIdx.x * blockDim.x + threadIdx.x;
int radius = FSIZE/2;
// 1.先初始化padding数据
//前raduis个线程分别初始化raduis长度的padding中的一个元素
if(threadIdx.x < radius)
{
//判断当前块是否是一个块
if(blockIdx.x == 0)
shareM[threadIdx.x] = 0;
else
//这里使用块的临近padding的前raduis个线程做padding初始化；也可以在非最后一个块中用前raduis个线程做右侧padding的初始化，
//写法为 shareM[blockDim.x + radius + threadIdx.x] = in[(blockIdx.x+1) * blockDim.x + threadIdx]
shareM[threadIdx.x] = in[id-radius];
}
else if(threadIdx.x >= blockDim.x-radius && threadIdx.x < blockDim.x)
{
//判断当前是否最后一个块
if(blockIdx.x == (len + blockDim.x - 1) / blockDim.x -1)
shareM[threadIdx.x + 2*radius] = 0;
else
shareM[threadIdx.x + 2*radius] = in[id+radius];
}
// 2. 初始化非padding数据
//初始化本线程对应输入计算的数据到共享内存
shareM[threadIdx.x + radius] = in[id];
// 3. 线程同步
__syncthreads();
// 4. 计算
double result = 0.0;
for (int i = 0;i < FSIZE;i++)
{
int offPos = threadIdx.x + i;
if (offPos >= 0 && offPos < len)
result += fliterConst[i] * shareM[offPos];
}
out[id] = result;
}
//串行版本
float serial()
{
//host串行计算
long start = getTime();
serialConv(in, out, QW, filter, FSIZE);
long end = getTime();
return end-start;
}
//GPU处理
float gpu()
{
//定义gpu全局内存数据空间
double* inD=NULL;
double* outD=NULL;
double* filterD=NULL;
//这里必须要取地址 &inD, 否则开辟空间失败
printf("memory in err %d\n",cudaMalloc((void**)&inD, QW * sizeof(double)));
printf("memory out err %d\n",cudaMalloc((void**)&outD, QW * sizeof(double)));
printf("memory filter err %d\n",cudaMalloc((void**)&filterD, FSIZE* sizeof(double)));
//拷贝数据
cudaMemcpy(inD, in, QW * sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(filterD, filter, FSIZE * sizeof(double), cudaMemcpyHostToDevice);
//核函数计算
//定义blockDim
dim3 blockDim(BlockSize,1,1);
//grid可以使用一维
dim3 gridDim((QW + blockDim.x - 1) / blockDim.x, 1, 1);
//cuda 记录时间,单位为 ms
cudaEvent_t startD, endD;
float gpuTime = 0.0;
cudaEventCreate(&startD);
cudaEventCreate(&endD);
cudaEventRecord(startD, 0);
kernalConv<<<gridDim, blockDim>>>(inD,outD, QW, filterD, FSIZE);
cudaEventRecord(endD, 0);
//等待事件完成
cudaEventSynchronize(endD);
cudaEventElapsedTime(&gpuTime, startD, endD);
cudaEventDestroy(startD);
cudaEventDestroy(endD);
//拷贝计算结果
cudaMemcpy(outFD, outD, QW * sizeof(double), cudaMemcpyDeviceToHost);
cudaFree(inD);
cudaFree(outD);
cudaFree(filterD);
return gpuTime;
}
//对卷积核使用常量内存优化，
float const_opt()
{
//定义gpu全局内存数据空间
double* inD=NULL;
double* outD=NULL;
// double* filterD=NULL;
//这里必须要取地址 &inD, 否则开辟空间失败
printf("memory in err %d\n",cudaMalloc((void**)&inD, QW * sizeof(double)));
printf("memory out err %d\n",cudaMalloc((void**)&outD, QW * sizeof(double)));
// printf("memory filter err %d\n",cudaMalloc((void**)&filterD, FSIZE* sizeof(double)));
//拷贝数据
cudaMemcpy(inD, in, QW * sizeof(double), cudaMemcpyHostToDevice);
// cudaMemcpy(filterD, filter, FSIZE * sizeof(double), cudaMemcpyHostToDevice);
//常量数据拷贝
cudaMemcpyToSymbol(fliterConst,filter,FSIZE*sizeof(double));
//核函数计算
//定义blockDim
dim3 blockDim(BlockSize,1,1);
//grid可以使用一维
dim3 gridDim((QW + blockDim.x - 1) / blockDim.x, 1, 1);
//cuda 记录时间,单位为 mse
cudaEvent_t startD, endD;
float gpuTime = 0.0;
cudaEventCreate(&startD);
cudaEventCreate(&endD);
cudaEventRecord(startD, 0);
kernalConv_const<<<gridDim, blockDim>>>(inD,outD, QW);
cudaEventRecord(endD, 0);
//等待事件完成
cudaEventSynchronize(endD);
cudaEventElapsedTime(&gpuTime, startD, endD);
cudaEventDestroy(startD);
cudaEventDestroy(endD);
//拷贝计算结果
cudaMemcpy(outFD, outD, QW * sizeof(double), cudaMemcpyDeviceToHost);
cudaFree(inD);
cudaFree(outD);
return gpuTime;
}
float const_shared_opt()
{
//定义gpu全局内存数据空间
double* inD=NULL;
double* outD=NULL;
//这里必须要取地址 &inD, 否则开辟空间失败
printf("memory in err %d\n",cudaMalloc((void**)&inD, QW * sizeof(double)));
printf("memory out err %d\n",cudaMalloc((void**)&outD, QW * sizeof(double)));
//拷贝数据
cudaMemcpy(inD, in, QW * sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpyToSymbol(fliterConst, filter, FSIZE * sizeof(double));
//核函数计算
//定义blockDim
dim3 blockDim(BlockSize,1,1);
//grid可以使用一维
dim3 gridDim((QW + blockDim.x - 1) / blockDim.x, 1, 1);
//cuda 记录时间,单位为 ms
cudaEvent_t startD, endD;
float gpuTime = 0.0;
cudaEventCreate(&startD);
cudaEventCreate(&endD);
cudaEventRecord(startD, 0);
kernelConv_const_share<<<gridDim, blockDim>>>(inD,outD, QW);
cudaEventRecord(endD, 0);
//等待事件完成
cudaEventSynchronize(endD);
cudaEventElapsedTime(&gpuTime, startD, endD);
cudaEventDestroy(startD);
cudaEventDestroy(endD);
//拷贝计算结果
cudaMemcpy(outFD, outD, QW * sizeof(double), cudaMemcpyDeviceToHost);
cudaFree(inD);
cudaFree(outD);
return gpuTime;
}
int main()
{
float cpuTime = 0.0;
float gpuTime = 0.0;
//定义均值滤波卷积核
for (int i = 0;i < FSIZE;i++)
*(filter + i) = (double)1.0 / FSIZE;
//读取一维信号
readFile(in, QW);
//串行与GPU并行对比测试
cpuTime = serial();
//gpu并行
#if 1
gpuTime = gpu();
//计算结果对比
for (int i = 0; i< QW;i++)
{
if (fabs(out[i] - outFD[i]) > EP)
{
printf("host[%d] %2.16f. != device[%d] %2.16f\n", i, out[i], i, outFD[i]);
return 1;
}
}
memset(outFD,0,QW*sizeof(double));
printf("串行与gpu计算结果一致,cpu耗时 %f ms,gpu耗时 %f ms,加速比为 %f\n",cpuTime,gpuTime,cpuTime/gpuTime);
#endif
//对卷积核使用常量内存优化
#if 1
gpuTime = const_opt();
//计算结果对比
for (int i = 0; i< QW;i++)
{
if (fabs(out[i] - outFD[i]) > EP)
{
printf("host[%d] %2.16f. != device[%d] %2.16f\n", i, out[i], i, outFD[i]);
return 1;
}
}
memset(outFD,0,QW*sizeof(double));
printf("串行与const_opt计算结果一致,cpu耗时 %f ms,gpu耗时 %f ms,加速比为 %f\n",cpuTime,gpuTime,cpuTime/gpuTime);
#endif
//1. 卷积核使使用常量内存
//2. 参与block计算输入数据放在共享内存
#if 1
gpuTime = const_shared_opt();
//计算结果对比
for (int i = 0; i< QW;i++)
{
if (fabs(out[i] - outFD[i]) > EP)
{
printf("host[%d] %2.16f. != device[%d] %2.16f\n", i, out[i], i, outFD[i]);
return 1;
}
}
printf("串行与const_shared_opt计算结果一致,cpu耗时 %f ms,gpu耗时 %f ms,加速比为 %f\n",cpuTime,gpuTime,cpuTime/gpuTime);
#endif
//保存计算结果
saveFile(out, QW);
//释放空间
free(in);
free(out);
free(outFD);
free(filter);
return 0;
}
3. 测试结果
| 项目 | 耗时(ms) | 加速比 |
|---|---|---|
| 串行 | 215 | |
| gpu | 4.3 | 49.96 |
| gpu + 常量内存 | 1.96 | 109 |
| gpu + 常量内存 + 共享内存 | 1.54 | 139.6 |
4. 结果分析
延迟分析
(1) 使用常量内存能显著减少访存延迟;
(2) 使用共享内存也能极大的提升访存延迟。
访存次数分析
数据长度：len，卷积核长度: 2*radius+1 = FilterSize;
最初的base版本:
单个内部线程块的访存次数（包含卷积核和计算输入数据）为: blockDim.x * (2FilterSize);
单个边缘线程块的访存次数为: 原次数-单边padding减少次数( 项数(最多减少访问次数+最少减少访问次数)/2 * 2[卷积核和计算数据都减少，加倍] ):
blockDim.x * (2*FilterSize) - radius * (radius+1)/2 * 2;
卷积核使用常量内存存储:
单个内部线程块的访存次数为: blockDim.x * FilterSize;
单个边缘线程块的访存次数为: blockDim * FilterSize - raduis * (radius + 1)/2;
卷积核使用常量内存，输入计算数据使用共享内存时：
单个内部线程块的访存次数为: blockDim.x + 2 * radius;
单个边缘线程块的访存次数为: blockDim.x + radius
使用常量内存和常量内存+共享内存访存对比：
单个内部线程块访存对比：(blockDim.x * FilterSize)/(blockDim.x + 2 * radius)
单个边缘线程块的访存对比：(blockDim * FilterSize - raduis * (radius + 1)/2)/(blockDim.x + radius)
假设 blickDim.x远大于FilterSize,可以简化为：
单个内部线程块访存对比：(blockDim.x * FilterSize)/blockDim.x = FilterSize；
单个边缘线程块的访存对比：(blockDim * FilterSize)/blockDim.x = FilterSize；
浙公网安备 33010602011771号