高性能計算-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號