高性能计算-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;
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· DeepSeek R1 简明指南:架构、训练、本地部署及硬件要求
· NetPad:一个.NET开源、跨平台的C#编辑器
· PowerShell开发游戏 · 打蜜蜂