На этом шаге мы рассмотрим использование разделяемой памяти.

    Рассмотренный в предыдущем шаге вариант упирается в доступ к глобальной памяти - для каждого тела нужно произвести N чтений положений других тел. При этом если рассмотреть все нити одного блока, то видно, что они все обращаются к одним и тем же данным, то есть очень много избыточных чтений.

    Воспользуемя довольно стандартным для CUDA приемом - дадим каждому блоку массив (размером, равным размеру блока) в разделяемой памяти. Тогда работу ядра можно будет организовать следующим образом: весь исходный массив делится на тайлы (длиной, равной размеру блока). Далее ядро для каждого тайла загружает его в разделяемую память, добавляет влияние всех тел данного тайла на тела, просчитываемые нитями блока, после чего переходит к следующему тайлу.

    Поскольку размер массива в разделяемой памяти равен размеру блока, то чтобы его загрузить, достаточно, чтобы каждая нить загрузила всего по одному элементу. Поскольку нити не выполняются физически параллельно, то после того, как нить загрузит свой элемент, необходима синхронизация, чтобы убедиться, что все остальные нити блока также загрузили свои элементы. Аналогично, после того как нить добавит силы, вызванные телами текущего тайла, также нужна синхронизация перед загрузкой следующего тайла.

    Приведем ниже фрагмент кода, использующего разделяемую память.

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

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

 __shared__ float4 sp[BLOCK_SIZE]; //Выделяем место в разделяемой памяти

 for(int i = 0; i < N/BLOCK_SIZE; ++i, ind += BLOCK_SIZE){
  //Копируем в разделяемую память все элементы, соответствующие данному блоку
  sp[threadIdx.x] = oldPos[ind + threadIdx.x]; 
  //Дожидаемся, когда все элементы будут загружены
  __syncthreads();

  //Производим расчет полной силы
  for(int j = 0; j < 256; ++j){
   float3 r;
  
   r.x = sp[j].x - pos.x;
   r.y = sp[j].y - pos.y;
   r.z = sp[j].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;
  }

  //Дожидаемся, когда будет вычислена полная сила для всех элементов данного блока
  __syncthreads();
 }

 float4 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;
}

    Весь проект можно взять здесь.

    Приведем ниже результаты выполнения приложения при использовании различных типов памяти.

Таблица 1. Результаты выполнения приложения

Тип памяти	Время выполнения (в миллисекундах)
float3	3383.6032715
float4	3738.6459961
foat4 (разделяемая память)	3198.6081543

    На следующем шаге мы рассмотрим библиотеку Thrust.

Предыдущий шаг Содержание Следующий шаг