c - 在CUDA中并行处理for循环(1D天真卷积)

Question

有人可以帮我将嵌套的for循环转换为CUDA内核吗？
这是我尝试转换为CUDA内核的函数：

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

Answer 1

您的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循环迭代时，将以倾斜的方式连续计算上图中的每一行，每行一个。第二个内核需要（2 * N）-1个线程，每个列/ out值一个。

我的第二个实现（conv_Kernel2）免除了对 temp数组的需要，只为每个column / 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)); // 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-memcheck运行任何CUDA代码，以查看是否发生了API错误。 。

编辑：我尝试了我的 conv_Kernel2的一个简单的共享内存版本，但没有更快。我相信这样做的原因是这些数据集（在 N = 4096， A和 B分别为16Kbytes， out大约为32Kbytes）足够小，可以轻松地放入GPU L2缓存中，没有rash动。

但是，对于较新的体系结构（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
$

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


在命令行上指定 N
仅包含cpu conv和gpu conv_Kernel2。
GPU时序测量中包括将数据移至GPU或从GPU移出的时间成本
提供了 typedef ... mytype;，以便可以轻松地重新编译代码以测试各种数据类型的行为。
打印出“加速因子”，即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;
}