const int N = 256*512; //Количество тел
const int BLOCK_SIZE = 256; //Размер блока
const float EPS = 0.0001f; //Погрешность

/*
 Ядро - функция принимает 5 параметров: 
	1. массив, в котором содержатся новые позиции тел
	2. массив, в котором содержатся новые скорости тел
	3. массив, в котором содержатся старые позиции тел
	4. массив, в котором содержатся старые скорости тел
	5. промежуток времени
*/
__global__ void kernelFloat3(float3 *newPos, float3 * newVel, 
          float3 *oldPos, float3 *oldVel, float dt){
 int index = blockIdx.x*blockDim.x + threadIdx.x; //Определяем номер тела

 float3 pos = oldPos[index]; //Запоминаем старую позицию тела
 float3 f = make_float3(0.0f, 0.0f, 0.0f); 

 for(int i = 0; i < N; ++i){
  float3 pi = oldPos[i];
  //Определяем координаты вектора изменения положения
  float3 r = make_float3(pi.x - pos.x, pi.y - pos.y, pi.z - pos.z); 

  //Определяем расстояние между телами
  float invDist = 1.0f/sqrtf(r.x*r.x + r.y*r.y + r.z*r.z + EPS*EPS); 
  float s = invDist*invDist*invDist;

  //Изменяем каждую из координат вектора полной силы
  f.x += r.x*s; 
  f.y += r.y*s;
  f.z += r.z*s;
 }

 float3 vel = oldVel[index];
 //Изменяем каждую из координат вектора скорости тела
 vel.x += f.x*dt;
 vel.y += f.y*dt;
 vel.z += f.z*dt;

 //Изменяем координаты положения тела
 pos.x += vel.x*dt;
 pos.y += vel.y*dt;
 pos.z += vel.z*dt;

 //Запоминаем новую скорость и новое положение тела
 newPos[index] = pos;
 newVel[index] = vel;
}

//Функция для решения поставленной задачи с использованием типа float3
float solveFloat3(){
 float3 *p = new float3[N]; //Выделяем место в памяти под массив позиций
 float3 *v = new float3[N]; //Выделяем место в памяти под массив скоростей
 float3 *pDev[2] = {NULL, NULL};
 float3 *vDev[2] = {NULL, NULL};

 //Заполняем массив позиций и массив скоростей случайными данными
 for(int i = 0; i < N; ++i){
  p[i].x = rand()/(float)RAND_MAX - 0.5f;
  p[i].y = rand()/(float)RAND_MAX - 0.5f;
  p[i].z = rand()/(float)RAND_MAX - 0.5f;

  v[i].x = rand()/(float)RAND_MAX - 0.5f;
  v[i].y = rand()/(float)RAND_MAX - 0.5f;
  v[i].z = rand()/(float)RAND_MAX - 0.5f;
 }

 //Выделяем место в памяти видеокарты под массивы позиций и скоростей
 cudaMalloc((void **)&pDev[0], N*sizeof(float3));
 cudaMalloc((void **)&pDev[1], N*sizeof(float3));
 cudaMalloc((void **)&vDev[0], N*sizeof(float3));
 cudaMalloc((void **)&vDev[1], N*sizeof(float3));

 //Копируем из памяти компьютера в память видеокарты положения тел и их скорости
 cudaMemcpy(pDev[0], p, N*sizeof(float3), cudaMemcpyHostToDevice);
 cudaMemcpy(vDev[0], v, N*sizeof(float3), cudaMemcpyHostToDevice);

 //Определяем переменные для замера времени расчетов
 cudaEvent_t start, stop;
 float gpuTime = 0.0f;

 cudaEventCreate(&start);   //Создаем точку старта отсчета
 cudaEventCreate(&stop);    //Создаем точку завершения отсчета
 cudaEventRecord(start, 0); //Определяем, сколько времени прошло с момента 
                            //запуска программы

 //Запускаем функцию для расчетов
 int index = 0;
 for(int i = 0; i < 2; ++i, index ^= 1){
  kernelFloat3<<<dim3(N/BLOCK_SIZE), dim3(BLOCK_SIZE)>>>
                 (pDev[index^1], vDev[index^1], pDev[index], vDev[index], 0.01f);
 }
 
 //Определяем, сколько времени прошло с момента запуска программы
 cudaEventRecord(stop, 0); 
 //Убеждаемся, что все вычисления были закончены и синхронизируемся по точке отсчета
 cudaEventSynchronize(stop); 
 //В переменную gpuTime записываем разницу между точками start и stop
 cudaEventElapsedTime(&gpuTime, start, stop); 

 //Копируем рузельтаты вычисление из памяти видеокарты в память компьютера
 cudaMemcpy(p, pDev[index^1], N*sizeof(float3), cudaMemcpyDeviceToHost);
 cudaMemcpy(v, vDev[index^1], N*sizeof(float3), cudaMemcpyDeviceToHost);

 //Очищаем выделенную память на видеокарте
 cudaFree(pDev[0]);
 cudaFree(pDev[1]);
 cudaFree(vDev[0]);
 cudaFree(vDev[1]);

 cudaEventDestroy(start);
 cudaEventDestroy(stop);
 cudaDeviceReset();

 //Очищаем выделенную память на компьютере
 delete[] p;
 delete[] v;
 return gpuTime;
}