Потоки в CUDA могут выполняться в любом порядке. Вычисление ema[i-1], возможно, не было начато к тому моменту, когда вы пытаетесь вычислить ema[i] в другом потоке (это зависит от завершения вычисления 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. Если это преобразование кажется вам неясным, вы можете создать таблицу из первых нескольких записей последовательности, аналогично методике, описанной в в этом ответе .
Это работает, но каждый поток выполняет итерацию по всему входному массиву вплоть до его индекса. Будет один потокоблок (с самыми высокими индексами массива), который займет больше времени, чем все остальные. В худшем случае поток выполняет примерно ту же работу, что и последовательная версия, поэтому не очень интересная параллельная реализация.
Чтобы решить эту проблему, мы можем сделать замечание о нерекурсивном уравнении формы. Согласно вашему коду, термин (1.0 - k) всегда меньше 1, поскольку k равен 2, деленному на некоторое положительное целое число больше 2 (т.е. мы будем считать, что period равно 2 или больше). Поэтому термин (1.0 - k)^x со временем становится исчезающе малым по мере продолжения суммирования. Мы также предполагаем, что ваши данные ограничены в диапазоне, примерно как вы показали. В этом случае, когда суммирование продолжается, в конечном итоге суммируемые термины не оказывают заметного влияния на сумму суммы float. С этими допущениями мы свернем циклическую обработку, когда наш (1.0 - k)^x член станет достаточно маленьким, чтобы не оказать существенного влияния на результат.
С этими допущениями и модификациями мы можем создать код 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
$
Tesla V100, CUDA 10.1
Повторим, действительность вышеприведенного кода с повышением производительности при досрочном выходе зависит от входных данных с ограниченным диапазоном. Я не буду пытаться тщательно охватить статистику, но если вы не имеете представления о статистике ваших входных данных, то приведенный выше метод может оказаться недействительным.