在CUDA上，向量步进加法速度较慢

Question

我正在尝试在CUDA C++代码中运行向量步骤加法函数，但对于尺寸为5000000的大型浮点数组，它比我的CPU版本运行得更慢。以下是相关的CUDA和CPU代码：

#define THREADS_PER_BLOCK 1024
typedef float real;
__global__ void vectorStepAddKernel2(real*x, real*y, real*z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i < size)
    {
        x[i*xstep] = alpha* y[i*ystep] + beta*z[i*zstep];
    }
}

cudaError_t vectorStepAdd2(real *x, real*y, real* z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{

    cudaError_t cudaStatus;
    int threadsPerBlock = THREADS_PER_BLOCK;
    int blocksPerGrid = (size + threadsPerBlock -1)/threadsPerBlock;
    vectorStepAddKernel2<<<blocksPerGrid, threadsPerBlock>>>(x, y, z, alpha, beta, size, xstep, ystep, zstep);

    // cudaDeviceSynchronize waits for the kernel to finish, and returns
    // any errors encountered during the launch.
    cudaStatus = cudaDeviceSynchronize();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching vectorStepAddKernel!\n", cudaStatus);
        exit(1);
    }

    return cudaStatus;
}

//CPU function:

void vectorStepAdd3(real *x, real*y, real* z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{
    for(int i=0;i<size;i++)
    {
        x[i*xstep] = alpha* y[i*ystep] + beta*z[i*zstep];
    }
}

当三个数组的大小为5000000且size=50000（即以这种逐步方式相加50,000个元素）时，调用vectorStepAdd2比调用vectorStepAdd3计算速度慢。

有什么办法可以加快GPU代码的运行速度吗？ 我的设备是Tesla M2090 GPU。

谢谢。

- assassin

---

跨步访问并不适合GPU的内存子系统，它更喜欢连续访问。如果跨度很小（例如<10个元素）且向量很长，则通过纹理访问只读数组可能会有帮助，值得一试。如果您要构建的是sm_35平台，可以对函数原型进行简单更改，使您的代码自动利用纹理路径，通过LDG指令：vectorStepAddKernel2(real * __restrict__ x, const real * __restrict__ y, const real * __restrict__ z, ...)。 - njuffa

你使用的xstep、ystep和zstep的值是多少？ - talonmies

@talonmies - 我正在使用的xstep、ystep、zstep值分别为4、5、7...但是，它们作为参数动态传递给函数（正如您所看到的），因此可以是任何值。 - assassin

@njuffa - 对不起，我没有完全理解你的意思。我在构建架构特定代码方面没有太多经验，但我认为我的Visual Studio设置已经设置为在Itanium平台上构建x64架构。我认为指令是基于处理器而不是芯片组的，请纠正我如果我错了。 - assassin

x64 是您主机系统的 64 位 x86 架构。sm_35 = 计算功能 3.5。当使用 nvcc 编译此平台时，这是您传递的内容：nvcc -arch=sm_35 -o [obj-file] [src-file]。LDG 指令在具有 sm_35 架构的 GPU 上是新的，例如 K20、K20x 和 Geforce Titan。 - njuffa

Answer 1

回答您的问题“我该怎么做才能加快GPU代码的速度？”

首先，让我说明一下，所提出的操作X = alpha * Y + beta * Z每字节数据传输所需的计算强度并不大。因此，在这个特定的代码上，我无法打败CPU时间。然而，可以介绍两个加速此代码的方法：

1. 使用页锁定内存进行数据传输操作。对于GPU版本，数据传输时间减少了约2倍，这在GPU版本中占据了整个执行时间的主导地位。

2. 使用带跨步的复制技术和cudaMemcpy2D，如@njuffa在这里所建议的。结果是双重的：我们可以将数据传输量减少到仅需要进行计算的量，并且我们可以按照评论中所建议的方式对数据进行连续操作。这使得数据传输时间再次减少了约3倍，计算时间提高了约10倍。

这段代码提供了这些操作的示例：

#include <stdio.h>
#include <stdlib.h>


#define THREADS_PER_BLOCK 1024
#define DSIZE 5000000
#define WSIZE 50000
#define XSTEP 47
#define YSTEP 43
#define ZSTEP 41
#define TOL 0.00001f


#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

typedef float real;

__global__ void vectorStepAddKernel2(real *x, real *y, real *z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i < size)
    {
        x[i*xstep] = alpha* y[i*ystep] + beta*z[i*zstep];
    }
}

__global__ void vectorStepAddKernel2i(real *x, real *y, real *z, real alpha, real beta, int size)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i < size)
    {
        x[i] = alpha* y[i] + beta*z[i];
    }
}

void vectorStepAdd2(real *x, real *y, real *z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{

    int threadsPerBlock = THREADS_PER_BLOCK;
    int blocksPerGrid = (size + threadsPerBlock -1)/threadsPerBlock;
    vectorStepAddKernel2<<<blocksPerGrid, threadsPerBlock>>>(x, y, z, alpha, beta, size, xstep, ystep, zstep);
    cudaDeviceSynchronize();
    cudaCheckErrors("kernel2 fail");
}


void vectorStepAdd2i(real *x, real *y, real *z, real alpha, real beta, int size)
{

    int threadsPerBlock = THREADS_PER_BLOCK;
    int blocksPerGrid = (size + threadsPerBlock -1)/threadsPerBlock;
    vectorStepAddKernel2i<<<blocksPerGrid, threadsPerBlock>>>(x, y, z, alpha, beta, size);
    cudaDeviceSynchronize();
    cudaCheckErrors("kernel3 fail");
}

//CPU function:

void vectorStepAdd3(real *x, real*y, real* z, real alpha, real beta, int size, int xstep, int ystep, int zstep)
{
    for(int i=0;i<size;i++)
    {
        x[i*xstep] = alpha* y[i*ystep] + beta*z[i*zstep];
    }
}

int main() {

  real *h_x, *h_y, *h_z, *c_x, *h_x1;
  real *d_x, *d_y, *d_z, *d_x1, *d_y1, *d_z1;

  int dsize = DSIZE;
  int wsize = WSIZE;
  int xstep = XSTEP;
  int ystep = YSTEP;
  int zstep = ZSTEP;
  real alpha = 0.5f;
  real beta = 0.5f;
  float et;

/*
  h_x = (real *)malloc(dsize*sizeof(real));
  if (h_x == 0){printf("malloc1 fail\n"); return 1;}
  h_y = (real *)malloc(dsize*sizeof(real));
  if (h_y == 0){printf("malloc2 fail\n"); return 1;}
  h_z = (real *)malloc(dsize*sizeof(real));
  if (h_z == 0){printf("malloc3 fail\n"); return 1;}
  c_x = (real *)malloc(dsize*sizeof(real));
  if (c_x == 0){printf("malloc4 fail\n"); return 1;}
  h_x1 = (real *)malloc(dsize*sizeof(real));
  if (h_x1 == 0){printf("malloc1 fail\n"); return 1;}
*/

  cudaHostAlloc((void **)&h_x, dsize*sizeof(real), cudaHostAllocDefault);
  cudaCheckErrors("cuda Host Alloc 1 fail");
  cudaHostAlloc((void **)&h_y, dsize*sizeof(real), cudaHostAllocDefault);
  cudaCheckErrors("cuda Host Alloc 2 fail");
  cudaHostAlloc((void **)&h_z, dsize*sizeof(real), cudaHostAllocDefault);
  cudaCheckErrors("cuda Host Alloc 3 fail");
  cudaHostAlloc((void **)&c_x, dsize*sizeof(real), cudaHostAllocDefault);
  cudaCheckErrors("cuda Host Alloc 4 fail");
  cudaHostAlloc((void **)&h_x1, dsize*sizeof(real), cudaHostAllocDefault);
  cudaCheckErrors("cuda Host Alloc 5 fail");


  cudaMalloc((void **)&d_x, dsize*sizeof(real));
  cudaCheckErrors("cuda malloc1 fail");
  cudaMalloc((void **)&d_y, dsize*sizeof(real));
  cudaCheckErrors("cuda malloc2 fail");
  cudaMalloc((void **)&d_z, dsize*sizeof(real));
  cudaCheckErrors("cuda malloc3 fail");
  cudaMalloc((void **)&d_x1, wsize*sizeof(real));
  cudaCheckErrors("cuda malloc4 fail");
  cudaMalloc((void **)&d_y1, wsize*sizeof(real));
  cudaCheckErrors("cuda malloc5 fail");
  cudaMalloc((void **)&d_z1, wsize*sizeof(real));
  cudaCheckErrors("cuda malloc6 fail");

  for (int i=0; i< dsize; i++){
    h_x[i] = 0.0f;
    h_x1[i] = 0.0f;
    c_x[i] = 0.0f;
    h_y[i] = (real)(rand()/(real)RAND_MAX);
    h_z[i] = (real)(rand()/(real)RAND_MAX);
    }


  cudaEvent_t t_start, t_stop, k_start, k_stop;
  cudaEventCreate(&t_start);
  cudaEventCreate(&t_stop);
  cudaEventCreate(&k_start);
  cudaEventCreate(&k_stop);
  cudaCheckErrors("event fail");

  // first test original GPU version

  cudaEventRecord(t_start);
  cudaMemcpy(d_x, h_x, dsize * sizeof(real), cudaMemcpyHostToDevice);
  cudaCheckErrors("cuda memcpy 1 fail");
  cudaMemcpy(d_y, h_y, dsize * sizeof(real), cudaMemcpyHostToDevice);
  cudaCheckErrors("cuda memcpy 2 fail");
  cudaMemcpy(d_z, h_z, dsize * sizeof(real), cudaMemcpyHostToDevice);
  cudaCheckErrors("cuda memcpy 3 fail");


  cudaEventRecord(k_start);
  vectorStepAdd2(d_x, d_y, d_z, alpha, beta, wsize, xstep, ystep, zstep);
  cudaEventRecord(k_stop);

  cudaMemcpy(h_x, d_x, dsize * sizeof(real), cudaMemcpyDeviceToHost);
  cudaCheckErrors("cuda memcpy 4 fail");
  cudaEventRecord(t_stop);
  cudaEventSynchronize(t_stop);
  cudaEventElapsedTime(&et, t_start, t_stop);
  printf("GPU original version total elapsed time is: %f ms.\n", et);
  cudaEventElapsedTime(&et, k_start, k_stop);
  printf("GPU original kernel elapsed time is: %f ms.\n", et);

  //now test CPU version

  cudaEventRecord(t_start);
  vectorStepAdd3(c_x, h_y, h_z, alpha, beta, wsize, xstep, ystep, zstep);
  cudaEventRecord(t_stop);
  cudaEventSynchronize(t_stop);
  cudaEventElapsedTime(&et, t_start, t_stop);
  printf("CPU version total elapsed time is: %f ms.\n", et);
  for (int i = 0; i< dsize; i++)
    if (fabsf((float)(h_x[i]-c_x[i])) > TOL) {
      printf("cpu/gpu results mismatch at i = %d, cpu = %f, gpu = %f\n", i, c_x[i], h_x[i]);
      return 1;
      }


  // now test improved GPU version

  cudaEventRecord(t_start);
//  cudaMemcpy2D(d_x1, sizeof(real),  h_x, xstep * sizeof(real), sizeof(real), wsize, cudaMemcpyHostToDevice);
//  cudaCheckErrors("cuda memcpy 5 fail");
  cudaMemcpy2D(d_y1, sizeof(real),  h_y, ystep * sizeof(real), sizeof(real), wsize, cudaMemcpyHostToDevice);
  cudaCheckErrors("cuda memcpy 6 fail");
  cudaMemcpy2D(d_z1, sizeof(real),  h_z, zstep * sizeof(real), sizeof(real), wsize, cudaMemcpyHostToDevice);
  cudaCheckErrors("cuda memcpy 7 fail");

  cudaEventRecord(k_start);
  vectorStepAdd2i(d_x1, d_y1, d_z1, alpha, beta, wsize);
  cudaEventRecord(k_stop);

  cudaMemcpy2D(h_x1, xstep*sizeof(real), d_x1, sizeof(real), sizeof(real), wsize, cudaMemcpyDeviceToHost);
  cudaCheckErrors("cuda memcpy 8 fail");
  cudaEventRecord(t_stop);
  cudaEventSynchronize(t_stop);
  cudaEventElapsedTime(&et, t_start, t_stop);
  printf("GPU improved version total elapsed time is: %f ms.\n", et);
  cudaEventElapsedTime(&et, k_start, k_stop);
  printf("GPU improved kernel elapsed time is: %f ms.\n", et);

  for (int i = 0; i< dsize; i++)
    if (fabsf((float)(h_x[i]-h_x1[i])) > TOL) {
      printf("gpu/gpu improved results mismatch at i = %d, gpu = %f, gpu imp = %f\n", i, h_x[i], h_x1[i]);
      return 1;
      }

  printf("Results:i   CPU     GPU     GPUi \n");
  for (int i = 0; i< 20*xstep; i+=xstep)
    printf("    %d         %f      %f     %f    %f    %f\n",i, c_x[i], h_x[i], h_x1[i]);


  return 0;
}

如前所述，我仍然无法打败CPU时间，这要么是由于我自己缺乏编码技能，要么是因为此操作在GPU上基本上没有足够的计算复杂度，从而不具有足够的吸引力。尽管如此，以下是一些样本结果：

GPU original version total elapsed time is: 13.352256 ms.
GPU original kernel elapsed time is: 0.195808 ms.
CPU version total elapsed time is: 2.599584 ms.
GPU improved version total elapsed time is: 4.228288 ms.
GPU improved kernel elapsed time is: 0.027392 ms.
Results:i   CPU     GPU     GPUi
    0         0.617285      0.617285     0.617285
    47         0.554522      0.554522     0.554522
    94         0.104245      0.104245     0.104245
....

我们可以看到，改进后的内核与原始内核相比，总体上减少了约3倍，其中几乎全部是由于数据复制时间的减少。这种数据复制时间的减少是由于改进的2D memcpy使我们只需要复制实际使用的数据。(如果没有页面锁定内存，这些数据传输时间将增加一倍，大约为原来的两倍)。我们还可以看到，对于原始内核，内核计算时间大约比CPU计算快10倍，而对于改进后的内核，内核计算时间大约比CPU计算快100倍。然而，在考虑数据传输时间时，我们无法克服CPU速度。

最后一点要注意的是，cudaMemcpy2D操作的“成本”仍然很高。对于向量大小减小100倍，我们只看到时间减少了3倍。因此，跨步访问仍然是一种相对昂贵的GPU使用方式。如果我们只是传输50,000个连续元素的向量，我们预计将几乎线性地减少100倍的复制时间(与复制5000000个元素的向量相比)。这意味着复制时间将少于1毫秒，我们的GPU版本将比CPU更快，至少对于这个简单的单线程CPU代码而言。

- Robert Crovella

---

谢谢！实际上，我在我的机器上使用您改进的内核时观察到内核运行时间比CPU函数调用时间慢。我知道这可能是由于各种原因...但正如您指出的那样，似乎这不是一个足够具有挑战性的问题可以在GPU上解决（或者现代CPU非常快：）） - assassin