Динамически расширяющийся массив в ядре cuda

Question

Я пытаюсь запустить алгоритм Брандеса (в основном bfs с некоторыми дополнительными операциями и структурами данных) на GPU, и я выделяю каждому потоку вершину для запуска брандов. Проблема, с которой я сталкиваюсь, заключается в том, что в моем коде

мне нужно хранить родителей каждой вершины, которую посетили во время bfs

. В реализации CPU это очень легко сделать, создав карту вектора и вызывая push_back всякий раз, когда я нахожу нового родителя, который является технически динамически расширяемым массивом. Я понятия не имею, как это сделать в CUDA.

Вот пример кода для нужной мне функциональности:

    vector<int> distance;               //Initialized to 0
    vector<int> paths;                  //Initialized to 0
    vector<bool> visited;               //Initialized to false
    map <int, vector<int> > parents;    //Parent vector of each key is empty
    queue<int> q;


    // Running bfs from vertex
    q.push(vertex);                     
    while(!q.empty())
    {
        int source = q.front();
        q.pop();

        for(auto neighbour : adjacency_list[source])
        {
            if(!visited[neighbour])
            {
                visited[neighbour] = true;
                q.push(neighbour);
                distance[neighbour] = distance[source] + 1;
            }
            if(distance[neighbour] == distance[source] + 1)
            {
                paths[neighbour] += paths[source];
                parents[neighbour].push_back(source);
            }
        }
    }

    {
        // Use data accumulated above for calculations
        ....
    }

Это строка (функциональность), которую мне сложно реализовать в коде устройства

родители [сосед] .push_back (источник);

Мои впечатления:

1. Я мог бы перераспределить (максимальная степень графика) список родителей для каждой вершины, но это будет стоить мне много неиспользованной памяти

2. Хранить родительские отношения как ребра в массиве размером 2 * Края, но мне нужны все родители вершины вместе (хранящиеся непрерывно или в одном контейнере), что невозможно в этой реализации

3. Я знаю о куче памяти GPU, но не могу придумать, как использовать ее для моего использования

4. Наихудший сценарий: я сначала запускаю bfs, чтобы найти no. родителей для каждой вершины, а затем выделить соответствующую память для каждой, а затем снова запустить Брандес.

Answer 1

1. Я думаю, что ваше впечатление 1 может быть реализовано примерно так, как описано здесь (стек для каждого потока, предварительно выделенный). У него есть проблемы, которые вы упоминаете, связанные с перераспределением. В более новых графических процессорах память размером в несколько гигабайт (или больше) является обычной, поэтому проблема перераспределения может быть не очень серьезной, если общий объем памяти не является проблемой.

2. Я думаю, что ваше впечатление 2 может быть реализовано примерно так, как описано здесь (потокобезопасный векторный push-back для всего устройства). У него есть проблемы, о которых вы упомянули, связанные с отсутствием упорядочения результатов в векторе результатов. Возможно, их можно решить с помощью операции сортировки после завершения операции сбора.

(4. Похоже, у вас, вероятно, уже есть представление о том, как произвести «худшее» впечатление 4.)

3. Это оставляет впечатление 3. Мы могли бы использовать идею, которая представляет собой комбинацию впечатления 1 и впечатления 2, то есть создать вектор push_back для каждого потока, но использовать распределение по требованию через ядро ​​malloc или * 1018. *. Распределение памяти в ядре, как это довольно медленно, и не без собственных проблем (например, вам может потребоваться зарезервировать дополнительное пространство кучи, выделенная память кучи в ядре не может участвовать в передаче на хост, небольшие выделения могут быть неэффективными в памяти использование), но на самом деле нет никакого способа определить, какой из этих подходов может быть наилучшим, без дополнительной информации о размерах вашей проблемы. Если отслеживание родительских узлов является относительно редкой операцией при обходе графа, подход с динамическим распределением может не быть проблемой.

Вот пример того, как можно создать простой вектор (для потока):

$ cat t376.cu
#include <iostream>
#include <cstdio>

#include <assert.h>
template <typename T>
class cu_vec{  // simple implementation of per-thread "vector"
  const size_t alloc_block_size = 4096; // tuning parameter
  T *my_ptr;
  size_t n_items;
  size_t alloc_blocks;
  public:
    __host__ __device__
    cu_vec(){
      assert(sizeof(T) <= alloc_block_size);
      n_items = 0;
      my_ptr = (T *)new char[alloc_block_size];
      assert(my_ptr != NULL);
      alloc_blocks = 1;}

    __host__ __device__
    cu_vec(size_t sz){
      assert(sizeof(T) <= alloc_block_size);
      n_items = sz;
      alloc_blocks = (n_items*sizeof(T)+alloc_block_size-1)/alloc_block_size;
      my_ptr = (T *)new char[alloc_blocks*alloc_block_size];
      assert(my_ptr != NULL);
      memset(my_ptr, 0, alloc_blocks*alloc_block_size);}

    __host__ __device__
    ~cu_vec(){
      if (my_ptr != NULL) delete[] my_ptr;
      }

    __host__ __device__
    void push_back(T const &item){ // first test if we can just store new item
      if ((n_items+1)*sizeof(T) > alloc_blocks*alloc_block_size){
        T *temp = (T *)new char[(alloc_blocks+1)*alloc_block_size];
        assert(temp != NULL);
        memcpy(temp, my_ptr, alloc_blocks*alloc_block_size);
        delete[] my_ptr;
        my_ptr = temp;
        alloc_blocks++;}
      my_ptr[n_items] = item;
      n_items++;}

    __host__ __device__
    size_t size(){
      return n_items;}

    __host__ __device__
    void clear(){
      n_items = 0;}

    __host__ __device__
    T& operator[](size_t idx){
      assert(idx < n_items);
      return my_ptr[idx];}

    __host__ __device__
    T& pop_back(){
      if (n_items > 0){
        n_items--;}
      return my_ptr[n_items];}

    __host__ __device__
    T* data(){
      return my_ptr;}

    __host__ __device__
    size_t storage_ratio(){
      return alloc_block_size/sizeof(T);}
};

struct ss
{
   unsigned x;
   float y;
};

__global__ void test(){

  cu_vec<ss> my_vec;
  ss temp = {threadIdx.x, 2.0f};
  my_vec.push_back(temp);
  assert(my_vec.size() == 1);
  assert(my_vec.storage_ratio() >= 1);
  ss temp2 = my_vec[0];
  printf("threadIdx.x: %u, ss.x: %u, ss.y: %f\n", threadIdx.x, temp2.x, temp2.y);
  temp.y = 3.0f;
  my_vec[0].x = temp.x;
  my_vec[0].y = temp.y;
  ss temp3 = my_vec.pop_back();
  printf("threadIdx.x: %u, ss.x: %u, ss.y: %f\n", threadIdx.x, temp3.x, temp3.y);
  my_vec.clear();
  temp.x = 0;
  for (int i = 0; i < 10000; i++){
    my_vec.push_back(temp);
    temp.x++;}
  temp.x--;
  for (int i = 0; i < 10000; i++) {
    assert(my_vec.pop_back().x == temp.x);
    temp.x--;}
  cu_vec<ss> my_vec2(2);
  assert(my_vec2[1].x == 0);
  assert(my_vec2[1].y == 0.0f);
}

int main(){

  //default heap space is 8MB, if needed reserve more with:
  cudaDeviceSetLimit(cudaLimitMallocHeapSize, (1048576*32));
  test<<<1, 4>>>();
  cudaDeviceSynchronize();
}
$ nvcc -std=c++11 -o t376 t376.cu
$ cuda-memcheck ./t376
========= CUDA-MEMCHECK
threadIdx.x: 0, ss.x: 0, ss.y: 2.000000
threadIdx.x: 1, ss.x: 1, ss.y: 2.000000
threadIdx.x: 2, ss.x: 2, ss.y: 2.000000
threadIdx.x: 3, ss.x: 3, ss.y: 2.000000
threadIdx.x: 0, ss.x: 0, ss.y: 3.000000
threadIdx.x: 1, ss.x: 1, ss.y: 3.000000
threadIdx.x: 2, ss.x: 2, ss.y: 3.000000
threadIdx.x: 3, ss.x: 3, ss.y: 3.000000
========= ERROR SUMMARY: 0 errors
$

Код не был проверен больше, чем вы видите здесь.