Странное поведение ядра CUDA, генерирование случайных значений

Я пытаюсь смоделировать систему пружинных масс с помощью CUDA. Ниже приведено ядро, которое обновляет положение частиц:

__global__ void timestep(float3 *pos, float3 *pos_antiga, float3 *aceleracao, int numParticulas) {

int index = threadIdx.x + blockIdx.x * blockDim.x;

if(index > (numParticulas - 1))
return;
t
float3 temp = pos[index];

pos[index].x = pos[index].x + (pos[index].x - pos_antiga[index].x) * (1.0f - DAMPING) + aceleracao[index].x * TIMESTEP;
pos[index].y = pos[index].y + (pos[index].y - pos_antiga[index].y) * (1.0f - DAMPING) + aceleracao[index].y * TIMESTEP;
pos[index].z = pos[index].z + (pos[index].z - pos_antiga[index].z) * (1.0f - DAMPING) + aceleracao[index].z * TIMESTEP;

pos_antiga[index] = temp;

}

pos представляет трехмерный вектор фактической позиции, pos_antiga позиция на предыдущем шаге по времени, DAMPING 0,01 и TIMESTEP 0,25. Я использую интеграцию Verlet. В системе без какой-либо силы, aceleracao ноль, поэтому pos а также pos_antigo одинаковы до и после вызовов ядра.

Однако после первой итерации CUDA возвращает сумасшедшие значения для некоторых осей, таких как 1.QNAN и -1.6241e + 016. Я думаю, что это связано с размерами блоков и потоков. Вызов ядра ниже:

timestep<<<16, 16>>>(pos_d, pos_antiga_d, aceleracao_d, numParticulas);

Итак, что мне не хватает?

РЕДАКТИРОВАТЬ: Ниже код вызывающего абонента:

void timestepGPU(vector<Particula> *p) {

// vector<Particula> has all the particles of the system.

// CPU
float *pos;
float *pos_antiga;
float *aceleracao;

// GPU
float *pos_d;
float *pos_antiga_d;
float *aceleracao_d;

// Number of particles
int numParticulas = p->size();

// Init
pos = new float[numParticulas * 3];
pos_antiga = new float[numParticulas * 3];
aceleracao = new float[numParticulas * 3];

// Transfering the values from the class to a plain vector
vector<Particula>::iterator p_tmp;
int i = 0;
for(p_tmp = p->begin(); p_tmp != p->end(); p_tmp++)
{
pos[i] = (*p_tmp).getPos().f[0];
pos[i + 1] = (*p_tmp).getPos().f[1];
pos[i + 2] = (*p_tmp).getPos().f[2];

pos_antiga[i] = (*p_tmp).getPosAntiga().f[0];
pos_antiga[i + 1] = (*p_tmp).getPosAntiga().f[1];
pos_antiga[i + 2] = (*p_tmp).getPosAntiga().f[2];

aceleracao[i] = (*p_tmp).getAceleracao().f[0];
aceleracao[i + 1] = (*p_tmp).getAceleracao().f[1];
aceleracao[i + 2] = (*p_tmp).getAceleracao().f[2];

i += 3;
}

// Here, I print the particle data BEFORE moving it to GPU
cout << "PRINT PARTICLE DATA" << endl;
for(i = 0; i < numParticulas * 3; i += 3) {
cout << i/3 << " - Pos: " << pos[i] << " " << pos[i + 1] << " " << pos[i + 2] << " | Pos Ant: " << pos_antiga[i] << " " << pos_antiga[i + 1] << " " << pos_antiga[i + 2] << " | Acel: " << aceleracao[i] << " " << aceleracao[i + 1] << " " << aceleracao[i + 2] << endl;
}
cout << "END" << endl;

// GPU
ErroCUDA(cudaMalloc((void**) &pos_d, numParticulas * 3 * sizeof(float)));
ErroCUDA(cudaMalloc((void**) &pos_antiga_d, numParticulas * 3 * sizeof(float)));
ErroCUDA(cudaMalloc((void**) &aceleracao_d, numParticulas * 3 * sizeof(float)));

// Moving data
ErroCUDA(cudaMemcpy(pos_d, pos, numParticulas * 3 * sizeof(float), cudaMemcpyHostToDevice));
ErroCUDA(cudaMemcpy(pos_antiga_d, pos_antiga, numParticulas * 3 * sizeof(float), cudaMemcpyHostToDevice));
ErroCUDA(cudaMemcpy(aceleracao_d, aceleracao, numParticulas * sizeof(float), cudaMemcpyHostToDevice));

// Setting number of blocks and threads per block
unsigned int numThreads, numBlocos;
calcularGrid(numParticulas, 64, numBlocos, numThreads);
//cout << numBlocos << "----------" << numThreads << endl;

// Kernel
timestep<<<numBlocos, numThreads>>>((float3 *) pos_d, (float3 *) pos_antiga_d, (float3 *) aceleracao_d, numParticulas);
ErroCUDA(cudaPeekAtLastError());
cudaDeviceSynchronize();

// Moving data back to the CPU
ErroCUDA(cudaMemcpy(pos, pos_d, numParticulas * 3 * sizeof(float), cudaMemcpyDeviceToHost));
ErroCUDA(cudaMemcpy(pos_antiga, pos_antiga_d, numParticulas * 3 * sizeof(float), cudaMemcpyDeviceToHost));

// Printing the particles' data AFTER Kernel call. At my GT 4xx, close to the 48th particle, it starts to show crazy values
cout << "PARTICLE DATA" << endl;
for(i = 0; i < numParticulas * 3; i += 3) {
cout << i/3 << " - Pos: " << pos[i] << " " << pos[i + 1] << " " << pos[i + 2] << " | Pos Ant: " << pos_antiga[i] << " " << pos_antiga[i + 1] << " " << pos_antiga[i + 2] << " | Acel: " << aceleracao[i] << " " << aceleracao[i + 1] << " " << aceleracao[i + 2] << endl;
}
cout << "END" << endl;

system("pause");

i = 0;
for(p_tmp = p->begin(); p_tmp != p->end(); p_tmp++)
{
if((*p_tmp).getMovel())
{
(*p_tmp).setPos(Vetor(pos[i], pos[i + 1], pos[i + 2]));
(*p_tmp).setPosAntiga(Vetor(pos_antiga[i], pos_antiga[i + 1], pos_antiga[i + 2]));
(*p_tmp).setAceleracao(Vetor(0, 0, 0));
}

i += 3;
}

ErroCUDA(cudaFree(pos_d));
ErroCUDA(cudaFree(pos_antiga_d));
ErroCUDA(cudaFree(aceleracao_d));

free(pos);
free(pos_antiga);
free(aceleracao);
}

В моем примере атрибут p имеет 100 элементов (10 x 10 частиц). Это сетка, начинающаяся с (0, 0, 0) и переходящая к (20, 20, 20) в трехмерном пространстве.

Еще раз спасибо всем за помощь!