CUDA C/C++ 访问前一个数组元素答案

【问题标题】：CUDA C/C++ accessing previous array elementCUDA C/C++ 访问前一个数组元素
【发布时间】：2019-11-16 12:35:44
【问题描述】：

在学习 CUDA 时，我正在研究一个计算移动平均线的小项目。虽然我的简单移动平均线 (SMA) 运行良好（虽然速度很慢且未优化），但我的指数移动平均线 (EMA) 总是导致错误的数字。

我发现问题在于 *(ema + i - 1) 始终为 0。同样的数组访问概念在测试 C++ 文件中运行良好，但在我的 CUDA 应用程序中却不行。我假设我只是不知道一些关于指针或 CUDA 的概念。


using namespace std;

// simple_ma not included

void __global__ exponential_ma(int n, int period, float *data, float *ema){
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if(i == 0){
        *ema = *data;
    }else if(i < n){

        float k = 2.0f/(period+1);
        *(ema + i) = *(data + i)*k + *(ema + i - 1) * (1.0f-k);
        // PROBLEM OCCURS ON THE LINE ABOVE, neither does ema[i-1] work
    }
}

int main(){
    /**
     * Function that computes a moving average on a vector
     */
    int N = 1<<5; // data size
    cout << "N = " << N << " bytes = " << N*sizeof(float) << endl;
    int period = 10; // moving average period

    // malloc'ed for stack usage instead of small heap size
    float *data = (float*)malloc(N*sizeof(float));
    float *sma = (float*)malloc(N*sizeof(float));
    float *ema = (float*)malloc(N*sizeof(float));

    float *d_data; // device pointer for data
    float *d_sma; // device pointer for simple moving average
    float *d_ema; // device pointer for exponential moving average

    // CUDA allocate memory for data, SMA, and EMA
    cudaMalloc(&d_data, N*sizeof(float));
    cudaMalloc(&d_sma, N*sizeof(float));
    cudaMalloc(&d_ema, N*sizeof(float));

    // initialize data
    srand(time(0));
    data[0] = rand() % 100 + 50;
    for(int i = 1; i < N; i++){
        data[i] = data[i-1] + rand() % 11 - 5;
    }

    // copy data from host to device
    cudaMemcpy(d_data, data, N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_sma, sma, N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_ema, ema, N*sizeof(float), cudaMemcpyHostToDevice);

    // call device function
    simple_ma<<<(N+255)/256, 256>>>(N, period, d_data, d_sma);
    exponential_ma<<<(N+255)/256, 256>>>(N, period, d_data, d_ema);

    cudaMemcpy(sma, d_sma, N*sizeof(float), cudaMemcpyDeviceToHost);
    cudaMemcpy(ema, d_ema, N*sizeof(float), cudaMemcpyDeviceToHost);

    for(int i = 0; i < N; i += 1){
        cout << "i = " << i << " data = "<< data[i] << " ---sma---> " << sma[i] << " ---ema---> " << ema[i] << endl;
    }

    cudaFree(d_data);
    cudaFree(d_sma);
    cudaFree(d_ema);

    return 0;
}

【问题讨论】：

cuda 中的线程可以按任意顺序执行。当您尝试在另一个线程中计算 ema[i] 时，ema[i-1] 的计算可能尚未开始（这取决于 ema[i-1] 的计算是否完成）。您在此处尝试计算的内容可以表示为a recurrence relation，并且需要特定的技术来解决。您用于算法的简单串行实现的方法不会以线程并行方式工作。
嗯，是的，这是有道理的。在这种情况下，我会假设某些值不会为 0（因为顺序是随机的，肯定会计算至少一个先前的元素？）
谢谢罗伯特。我必须想办法让这项工作以线程并行的方式工作。关于如何做这样的事情有什么建议吗？
我的建议在我给出的链接答案中。我不确定它会让你到达那里。
This paper 似乎给出了并行计算 EMA 的方法的粗略草图。但是，它似乎依赖于全局同步。通过一些搜索，您可能会找到更好的论文。

标签： c++ arrays pointers cuda malloc

【解决方案1】：

CUDA 中的线程可以按任何顺序执行。当您尝试在另一个线程中计算 ema[i] 时，ema[i-1] 的计算可能尚未开始（这取决于 ema[i-1] 的计算是否完成）。您用于此算法的简单串行实现的方法不会以线程并行方式工作

考虑到这一点，这是一种可能的方法。

首先，重铸你的递归 ema 计算：

     ema[0] = data[0]
i>0: ema[i] = k*data[i]+(1-k)*ema[i-1]

非递归形式：

     ema[0] = data[0]
i>0: ema[i] = ((1-k)^i)*data[0] + ∑(((1-k)^x)*k*data[i-x])  
                                  x=0..i-1

这将告诉我们如何编写我们的 CUDA 内核代码。如果您觉得这种转换不太清楚，您可能希望创建一个包含序列前几个条目的表格，类似于this answer 中描述的方法。

它可以工作，但是每个线程都在遍历整个输入数组直到它的索引。将有一个线程块（具有最高数组索引）比所有其他线程块花费更长的时间。最坏情况下的线程所做的工作与串行版本大致相同，因此不是一个非常有趣的并行实现。

为了解决这个问题，我们可以对非递归形式方程进行观察。根据您的代码，术语 (1.0 - k) 始终小于 1，因为 k 是 2 除以某个大于 2 的正整数（即我们假设 period 是 2 或更大）。因此，随着求和的进行，(1.0 - k)^x 项最终会变得非常小。我们还将假设您的数据是有界的，大致如您所示。在这种情况下，随着求和的进行，最终求和的项对float sum 数量没有明显影响。有了这些假设，当我们的(1.0 - k)^x 项变得足够小以至于对结果没有显着影响时，我们将减少循环处理。

通过这些假设和修改，我们可以创建一个运行速度比原始串行 CPU 版本更快的 CUDA 代码，同时保持较小的误差范围。

$ cat t1444.cu
#include <iostream>
#include <cstdio>
#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)
#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}


__global__ void gpu_ema(const int n, const float k, const float * __restrict__ data, float * __restrict__ ema, const float tol){
  int i = blockIdx.x*blockDim.x+threadIdx.x;
  if (i == 0) ema[0] = data[0];
  else if (i < n){
    float sum = 0;
    float fac = 1.0f - k;
    float m = 1.0f;
    int j;
    for (j = 0; j < i; j++){
      sum += m*k*data[i-j];
      m *= fac;
      if (m < tol) break;  // early exit validity depends on a variety of assumptions
      }
    if (j == i) sum += m*data[0];
    ema[i] = sum;
    }
}

void cpu_ema(int n, int period, float *data, float *ema){
  ema[0] = data[0];
  float k = 2.0f/(period+1);
  for (int i = 1; i < n; i++)
    ema[i] = data[i]*k + ema[i-1]*(1.0f-k);
}
int main(){
    /**
     * Function that computes a moving average on a vector
     */
    int N = 1<<20; // data size
    std::cout << "N = " << N << " bytes = " << N*sizeof(float) << std::endl;
    int period = 10; // moving average period

    // malloc'ed for stack usage instead of small heap size
    float *data = (float*)malloc(N*sizeof(float));
    float *ema = (float*)malloc(N*sizeof(float));
    float *gema = (float*)malloc(N*sizeof(float));

    float *d_data; // device pointer for data
    float *d_ema; // device pointer for exponential moving average

    // CUDA allocate memory for data, SMA, and EMA
    cudaMalloc(&d_data, N*sizeof(float));
    cudaMalloc(&d_ema, N*sizeof(float));

    // initialize data
    srand(time(0));
    data[0] = rand() % 100 + 50;
    for(int i = 1; i < N; i++){
        data[i] = data[i-1] + rand() % 11 - 5;
    }

    // copy data from host to device
    cudaMemcpy(d_data, data, N*sizeof(float), cudaMemcpyHostToDevice);

    // call device function
    long long gpu_t = dtime_usec(0);
    gpu_ema<<<(N+255)/256, 256>>>(N, 2.0f/(period+1), d_data, d_ema, 1e-7);
    cudaDeviceSynchronize();
    gpu_t = dtime_usec(gpu_t);
    long long cpu_t = dtime_usec(0);
    cpu_ema(N, period, data, ema);
    cpu_t = dtime_usec(cpu_t);
    if (N < 33)
      for (int i = 0; i < N; i++)
        std::cout << ema[i] << ",";
    std::cout << std::endl;

    cudaMemcpy(gema, d_ema, N*sizeof(float), cudaMemcpyDeviceToHost);
    cudaCheckErrors("some CUDA error");
    if (N < 33)
      for(int i = 0; i < N; i += 1)
          std::cout << gema[i] << ",";
    std::cout << std::endl;
    float max_err = fabs(gema[0] - ema[0])/ema[0];
    for (int i = 1; i < N; i++)
      max_err = max(max_err, fabs(gema[i] - ema[i])/ema[0]);
    std::cout << "max err: " << max_err*100.0 << "% final gpu: " << gema[N-1] << " final cpu: " << ema[N-1] << std::endl;
    std::cout << "cpu time: " << cpu_t/(float)USECPSEC << "s gpu time: " << gpu_t/(float)USECPSEC  << "s" << std::endl;
    cudaFree(d_data);
    cudaFree(d_ema);

    return 0;
}
$ nvcc -o t1444 t1444.cu
$ ./t1444
N = 1048576 bytes = 4194304


max err: 0.00218633% final gpu: 1311.38 final cpu: 1311.38
cpu time: 0.006346s gpu time: 0.000214s
$

特斯拉 V100、CUDA 10.1

重复一遍，上述代码的有效性与性能增强提前退出取决于范围有限的输入数据。我不会尝试仔细地涵盖统计信息，但是如果您不了解输入数据的统计信息，那么上述方法可能无效。

【讨论】：