高性能计算-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++)
{
// offpos 写在循环中的好处，生命周期短，减少线程寄存器占用数量
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数据；
// grid第一个块/最后一个块的共享内存前/后 padding raduis个0
// 其他块内第一个/最后一个线程padding 有效数据
// (1)如果是block索引为0 的块，左侧的padding应为0，否则应为前一个数据块的后radius个数据
// (2)如果是block索引为最后一个块，右侧的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[threadIdx.x + i];
}
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号