你的CPU
conv
函数似乎是这样做的(以
N
= 4为例):
A0B0 A0B1 A0B2 A0B3 + ^
A1B0 A1B1 A1B2 A1B3 + N
A2B0 A2B1 A2B2 A2B3 + rows
A3B0 A3B1 A3B2 A3B3 = v
out0 out1 out2 out3 out4 out5 out6
<- (2*N)-1 columns ->
你的卷积(对我来说)的特点在于它将两个长度相等的信号进行卷积。由于GPU喜欢处理“大”问题,这意味着N应该很大。然而,你的conv_Kernel实现中存在一个立即的问题,即它意味着块维度将用于索引A,线程维度将用于索引B。但是线程维度(threadIdx.x)对于当前CUDA GPU仅限于512或1024。这将使我们只能解决相当小的问题。
你的实现存在各种其他问题。其中一个问题是分配的共享内存大小不足以容纳i+j范围(可以从0->2*(N-1))。当然,这很容易修复,但更严重的问题是我没有看到一种方式将你的算术映射到任何类似于上面所述的期望模式。经过一段时间思考你的内核后,我放弃了它。
卷积问题有许多与之相关的研究,并且可以通过各种方式进行优化以适应GPU等大规模并行架构。因此,我将专注于两个非常简单的实现,这些实现基于上面的图表立即提出。
第一个实现的方法是重新创建上面的图表。我们将创建一个中间的
temp
数组来存储所有单个的AxBy乘积,在
conv_Kernel
中计算并存储这些乘积。然后,我们将启动第二个内核(
sum_Kernel
),它只是对
temp
数组的列进行求和,以产生各种
out
值。第一个内核需要
N
个线程,这些线程将按斜向方式迭代通过
N
个for循环迭代中的每一行计算每一行的上述图表。第二个内核需要(2*N)-1个线程,每个线程用于每个列/
out
值。
我的第二个实现方法(conv_Kernel2)不需要
temp
数组,而是为每个列/
out
值分配一个线程,并逐行遍历
N
行,逐行计算必要的乘积,并“即时”地对这些乘积进行求和。然后直接将求和结果存储在
out
数组中。
仅考虑计算,而不是数据移动/初始化所需的时间,在K20x GPU上,当
N
=512时,GPU实现开始比朴素的单线程CPU实现更快。第二个实现也值得称赞,因为唯一需要移动的数据是A、B和结果。第一个实现还需要分配并将
temp
数组初始化为全部零。
temp
数组的大小与
N
*
N
成正比,因此第二个实现的好处在于它不需要这个临时存储空间。
下面是一个完整的测试用例,运行并计时您提供的CPU实现以及我创建的两个略有不同的GPU实现:
$ cat t617.cu
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#define N 4096
#define RG 10
#define USECPSEC 1000000ULL
#define nTPB 256
void conv(int* A, int* B, int* out) {
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j)
out[i + j] += A[i] * B[j];
}
unsigned long long dtime_usec(unsigned long long prev){
timeval tv1;
gettimeofday(&tv1,0);
return ((tv1.tv_sec * USECPSEC)+tv1.tv_usec) - prev;
}
__global__ void conv_Kernel(int* A, int *B, int* temp) {
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < N){
int my_B = B[idx];
for (int i = 0; i < N; i++)
temp[idx + (i*2*N) + i] = my_B * A[i];
}
}
__global__ void sum_Kernel(int *temp, int *out){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < (2*N)-1){
int my_sum = 0;
for (int i = 0; i < N; i++) my_sum += temp[idx + (i*2*N)];
out[idx] = my_sum;}
}
__global__ void conv_Kernel2(int *A, int *B, int *out){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < (2*N)-1){
int my_sum = 0;
for (int i = 0; i < N; i++)
if (((idx < N) && (i <= idx)) || ((idx >= N) && (i > (idx-N)))) my_sum += A[i]*B[idx-i];
out[idx] = my_sum;
}
}
int main(){
int *h_A, *d_A, *h_result, *d_result, *result, *h_B, *d_B, *A, *B, *d_temp;
B = (int *)malloc(N*sizeof(int));
A = (int *)malloc(N*sizeof(int));
h_A = (int *)malloc(N*sizeof(int));
h_B = (int *)malloc(N*sizeof(int));
h_result = (int *)malloc(2*N*sizeof(int));
result = (int *)malloc(2*N*sizeof(int));
cudaMalloc(&d_B, N*sizeof(int));
cudaMalloc(&d_A, N*sizeof(int));
cudaMalloc(&d_result, 2*N*sizeof(int));
cudaMalloc(&d_temp, 2*N*N*sizeof(int));
for (int i=0; i < N; i++){
A[i] = rand()%RG;
B[i] = rand()%RG;
h_A[i] = A[i];
h_B[i] = B[i];}
for (int i=0; i < 2*N; i++){
result[i] = 0;
h_result[i] = 0;}
unsigned long long cpu_time = dtime_usec(0);
conv(A, B, result);
cpu_time = dtime_usec(cpu_time);
cudaMemcpy(d_A, h_A, N*sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, N*sizeof(int), cudaMemcpyHostToDevice);
cudaMemset(d_result, 0, 2*N*sizeof(int));
cudaMemset(d_temp, 0, 2*N*N*sizeof(int));
unsigned long long gpu_time = dtime_usec(0);
conv_Kernel<<<(N+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_temp);
sum_Kernel<<<((2*(N-1))+nTPB-1)/nTPB, nTPB>>>(d_temp, d_result);
cudaDeviceSynchronize();
gpu_time = dtime_usec(gpu_time);
cudaMemcpy(h_result, d_result, 2*N*sizeof(int), cudaMemcpyDeviceToHost);
for (int i = 0; i < 2*N; i++) if (result[i] != h_result[i]) {printf("mismatch at %d, cpu: %d, gpu %d\n", i, result[i], h_result[i]); return 1;}
printf("Finished. Results match. cpu time: %ldus, gpu time: %ldus\n", cpu_time, gpu_time);
cudaMemset(d_result, 0, 2*N*sizeof(int));
gpu_time = dtime_usec(0);
conv_Kernel2<<<((2*(N-1))+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_result);
cudaDeviceSynchronize();
gpu_time = dtime_usec(gpu_time);
cudaMemcpy(h_result, d_result, 2*N*sizeof(int), cudaMemcpyDeviceToHost);
for (int i = 0; i < 2*N; i++) if (result[i] != h_result[i]) {printf("mismatch2 at %d, cpu: %d, gpu %d\n", i, result[i], h_result[i]); return 1;}
printf("Finished. Results match. cpu time: %ldus, gpu2 time: %ldus\n", cpu_time, gpu_time);
return 0;
}
$ nvcc -arch=sm_35 -o t617 t617.cu
$ ./t617
Finished. Results match. cpu time: 69059us, gpu time: 3204us
Finished. Results match. cpu time: 69059us, gpu2 time: 1883us
$ nvcc -arch=sm_35 -O3 -o t617 t617.cu
$ ./t617
Finished. Results match. cpu time: 13750us, gpu time: 3214us
Finished. Results match. cpu time: 13750us, gpu2 time: 1886us
$
(请注意,即使只是使用-O3参数也会显着提高CPU代码执行效率)
正如我所提到的,我认为我的两个示例对于GPU代码也相当“幼稚”(例如,都没有使用共享内存),但它们可能会给您一些启示,让您开始着手。
为了简洁起见,我没有进行CUDA错误检查。 但是,我建议每当您在CUDA代码上遇到问题时,都要进行
适当的CUDA错误检查。 对于您的
conv_Kernel
,我认为它会指出一些错误(如果您尝试运行它)。 作为快速测试,您始终可以使用
cuda-memcheck
运行任何CUDA代码,以查看是否发生任何API错误。
编辑: 我尝试使用简单的共享内存版本来实现我的conv_Kernel2
,但是没有加快速度。我认为这是因为这些数据集(在N
=4096时,A
和B
每个16K字节,out
大约32K字节)足够小,可以轻松地适应GPU L2高速缓存,而不会出现抖动。
然而,对于新的架构(cc 3.5及更高版本),CUDA编译器有时可以将只读输入数据适当地标识为内核如果。因此,如果我们将我的conv_Kernel2
定义更改为:
__global__ void conv_Kernel2(const int * __restrict__ A, const int * __restrict__ B, int *out){
然后我看到执行时间略有改善,就我个人而言:
$ ./t617
Finished. Results match. cpu time: 13792us, gpu time: 3209us
Finished. Results match. cpu time: 13792us, gpu2 time: 1626us
$
我创建了一个修改版的代码,它执行以下操作:
- 在命令行上指定
N
- 只包括CPU的
conv
和GPU的conv_Kernel2
- 将数据从GPU传输到主机和反之所需的时间计入GPU计时测量中
- 提供了一个
typedef ... mytype;
,以便可以轻松重新编译代码以测试不同数据类型的行为。
- 输出“加速比”,即CPU时间除以GPU时间。
修改后的代码:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#define MAXN 100000
#define RG 10
#define USECPSEC 1000000ULL
#define nTPB 256
typedef double mytype;
void conv(const mytype *A, const mytype *B, mytype* out, int N) {
for (int i = 0; i < N; ++i)
for (int j = 0; j < N; ++j)
out[i + j] += A[i] * B[j];
}
unsigned long long dtime_usec(unsigned long long prev){
timeval tv1;
gettimeofday(&tv1,0);
return ((tv1.tv_sec * USECPSEC)+tv1.tv_usec) - prev;
}
__global__ void conv_Kernel2(const mytype * __restrict__ A, const mytype * __restrict__ B, mytype *out, const int N){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
if (idx < (2*N)-1){
mytype my_sum = 0;
for (int i = 0; i < N; i++)
if (((idx < N) && (i <= idx)) || ((idx >= N) && (i > (idx-N)))) my_sum += A[i]*B[idx-i];
out[idx] = my_sum;
}
}
int main(int argc, char *argv[]){
mytype *h_A, *d_A, *h_result, *d_result, *result, *h_B, *d_B, *A, *B;
if (argc != 2) {printf("must specify N on the command line\n"); return 1;}
int my_N = atoi(argv[1]);
if ((my_N < 1) || (my_N > MAXN)) {printf("N out of range\n"); return 1;}
B = (mytype *)malloc(my_N*sizeof(mytype));
A = (mytype *)malloc(my_N*sizeof(mytype));
h_A = (mytype *)malloc(my_N*sizeof(mytype));
h_B = (mytype *)malloc(my_N*sizeof(mytype));
h_result = (mytype *)malloc(2*my_N*sizeof(mytype));
result = (mytype *)malloc(2*my_N*sizeof(mytype));
cudaMalloc(&d_B, my_N*sizeof(mytype));
cudaMalloc(&d_A, my_N*sizeof(mytype));
cudaMalloc(&d_result, 2*my_N*sizeof(mytype));
for (int i=0; i < my_N; i++){
A[i] = rand()%RG;
B[i] = rand()%RG;
h_A[i] = A[i];
h_B[i] = B[i];}
for (int i=0; i < 2*my_N; i++){
result[i] = 0;
h_result[i] = 0;}
unsigned long long cpu_time = dtime_usec(0);
conv(A, B, result, my_N);
cpu_time = dtime_usec(cpu_time);
cudaMemset(d_result, 0, 2*my_N*sizeof(mytype));
unsigned long long gpu_time = dtime_usec(0);
cudaMemcpy(d_A, h_A, my_N*sizeof(mytype), cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, my_N*sizeof(mytype), cudaMemcpyHostToDevice);
conv_Kernel2<<<((2*(my_N-1))+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_result, my_N);
cudaDeviceSynchronize();
cudaMemcpy(h_result, d_result, 2*my_N*sizeof(mytype), cudaMemcpyDeviceToHost);
gpu_time = dtime_usec(gpu_time);
for (int i = 0; i < 2*my_N; i++) if (result[i] != h_result[i]) {printf("mismatch2 at %d, cpu: %d, gpu %d\n", i, result[i], h_result[i]); return 1;}
printf("Finished. Results match. cpu time: %ldus, gpu time: %ldus\n", cpu_time, gpu_time);
printf("cpu/gpu = %f\n", cpu_time/(float)gpu_time);
return 0;
}