2014-12-02 73 views
1

有人可以帮助我将嵌套for循环转换为CUDA内核吗? 这里是我想转换成CUDA内核的功能:在CUDA中并行化一个for循环(1D Naive Convolution)

// Convolution on Host 
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]; 
} 

我已经很努力并行的代码。
这里是我的尝试:

__global__ void conv_Kernel(int* A, int* B, int* out) { 

    int i = blockIdx.x; 
    int j = threadIdx.x; 

    __shared__ int temp[N]; 

    __syncthreads(); 
    temp[i + j] = A[i] * B[j]; 
    __syncthreads(); 

    int sum = 0; 
    for (int k = 0; k < N; k++) 
     sum += temp[k]; 
    out[i + j] = sum; 
} 

回答

2

你的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。但是对于当前的CUDA GPU,线程尺寸(threadIdx.x)限制为512或1024。这将使我们只能解决相当小的问题。

还有其他各种问题的认识。一个问题是分配的共享内存大小不足以适应i+j范围(可能从0-> 2 *(N-1))。这当然是微不足道的,但更严重的问题是,我没有看到一种方法将算法映射到类似上面所需的模式的任何东西上。花了一点时间思考你的内核后,我放弃了它。

卷积问题有与之相关联进行了大量研究,并可以在像GPU的大规模并行架构的各种方式进行优化。因此,我将重点介绍两个非常简单的实现,它们立即根据上图进行建议。

第一实现是简单地上图中重新创建。我们将创建一个中间temp阵列来存储所有单独的AxBy产品,计算并将这些产品存储在conv_Kernel中。然后,我们将启动第二个内核(sum_Kernel),它简单地求和temp数组的列,以产生各种值out值。第一个内核需要N线程,它将连续计算上面图的每一行,当我们迭代通过N for-loop迭代时,以斜线方式,每行一个。第二个内核需要(2 * N)-1个线程,每列一个/值。

我的第二个实现(conv_Kernel2)与需要一个temp阵列分配,和只有一个线程分配给每个列/ out值,并重复通过N行,计算所需的产品行接一行,和求和这些产品“即时”。总和结果然后直接存储在out阵列中。

考虑到只有计算,而不是数据移动/初始化所需的时间,GPU实现开始比K20x GPU上N = 512的天真单线程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)); // just for error checking, the kernel2 require no initialization of the result 

    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代码问题时,执行proper cuda error checking。对于您的conv_Kernel,我相信它会显示一些错误(如果您尝试运行它)。作为快速测试,您始终可以运行任何CUDA代码(cuda-memcheck)以查看是否发生任何API错误。

编辑:我试过我的conv_Kernel2简单的共享内存版本,但它没有更快。我相信其原因是这些数据集(在N = 4096,AB每个16KB,out约为32KB)足够小,可以轻松地安装到GPU L2缓存中,而不会发生抖动。但是,对于较新的体系结构(cc 3.5和更新的版本),CUDA编译器有时可以对内核进行额外的优化if the read-only input data is properly identified。因此,如果我们改变我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 
$ 

我创建了执行以下操作的代码的修改版本:

  1. N在命令行
  2. 只有CPU上指定包括和gpu conv_Kernel2
  3. 将数据从GPU转移到/时间成本包括在GPU定时测量
  4. 一个typedef ... mytype;被设置为使得可以将代码重新编译的容易来测试各种数据类型的行为。
  5. 打印出“加速因子”,即CPU时间除以GPU时间。

修改后的代码:

#include <stdio.h> 
#include <stdlib.h> 
#include <time.h> 
#include <sys/time.h> 

// RG*RG*MAXN must fit within mytype 

#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; 
}