Параллельные (GPU) алгоритмы для асинхронных клеточных автоматов

У меня есть коллекция вычислительных моделей, которые можно описать как асинхронные клеточные автоматы. Эти модели напоминают модель Изинга, но немного сложнее. Кажется, что такие модели выиграли бы от работы на GPU, а не на CPU. К сожалению, распараллелить такую ​​модель довольно непросто, и мне совершенно не понятно, как это сделать. Я знаю, что по этому вопросу есть литература, но все это, по-видимому, предназначено для компьютерных ученых, которые интересуются деталями алгоритмической сложности, а не для того, кто, как я, просто хочет описать что-то, что я могу реализовать, и следовательно, я нахожу это довольно непроницаемым.

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

Я также должен пояснить, что асинхронный клеточный автомат довольно сильно отличается от синхронного, и методы распараллеливания синхронных СА (такие как жизнь Конвея) не могут быть легко адаптированы к этой проблеме. Разница в том, что синхронный CA обновляет каждую ячейку одновременно на каждом временном шаге, тогда как асинхронный обновляет случайно выбранную локальную область на каждом временном шаге, как описано ниже.

Модели, которые я хочу распараллелить, реализованы на решетке (обычно гексагональной), состоящей из ~ 100000 ячеек (хотя я бы хотел использовать больше), и непараллельный алгоритм их запуска выглядит следующим образом:

1. Выберите соседнюю пару ячеек наугад

2. Вычислить «энергетическую» функцию на основе локальной окрестности, окружающей эти клетки$\Delta E$

3. С вероятностью, которая зависит от (с β параметром), либо поменяйте местами состояния двух ячеек, либо ничего не делайте.$e^{-\beta \Delta E}$$\beta$

4. Повторите вышеуказанные шаги до бесконечности.

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

Стоит отметить, что меня интересует не только состояние равновесия, но и динамика переходного процесса в этих системах, поэтому мне нужно нечто, имеющее динамику, эквивалентную приведенной выше, а не то, что будет приближаться к тому же равновесному распределению. (Таким образом, вариации алгоритма chequerboard не то, что я ищу.)

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

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

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

• Приблизьте систему следующим образом: пусть время протекает дискретно. Разделите решетку на другуюнабор областей на каждом временном шаге в соответствии с некоторой предварительно определенной схемой, и каждый блок графического процессора случайным образом выбирает и обновляет пару ячеек сетки, соседство которых не перекрывает границу региона. Поскольку границы меняются каждый раз по времени, это ограничение может не сильно влиять на динамику, если регионы относительно велики. Это кажется простым для реализации и, вероятно, будет быстрым, но я не знаю, насколько хорошо он будет приближаться к динамике, или какова лучшая схема выбора границ региона на каждом временном шаге. Я нашел некоторые ссылки на «блочно-синхронные клеточные автоматы», которые могут совпадать или не совпадать с этой идеей. (Я не знаю, потому что кажется, что все описания метода либо на русском языке, либо в источниках, к которым у меня нет доступа.)

Мои конкретные вопросы следующие:

• Являются ли какие-либо из вышеперечисленных алгоритмов разумным способом приблизиться к распараллеливанию GPU асинхронной модели CA?

• Есть ли способ лучше?

• Существует ли существующий библиотечный код для этого типа проблемы?

• Где я могу найти четкое англоязычное описание метода "block-синхронный"?

Прогресс

Я полагаю, что придумала способ распараллеливания асинхронного ЦС, который мог бы подойти. Алгоритм, описанный ниже, предназначен для обычного асинхронного ЦС, который обновляет только одну ячейку за раз, а не для соседней пары ячеек, как моя. Есть некоторые проблемы с обобщением этого в моем конкретном случае, но я думаю, что у меня есть идея, как их решить. Тем не менее, я не уверен, какую выгоду от скорости это даст по причинам, обсуждаемым ниже.

Идея состоит в том, чтобы заменить асинхронный CA (далее ACA) стохастическим синхронным CA (SCA), который ведет себя эквивалентно. Для этого сначала представим, что ACA - это пуассоновский процесс. То есть время протекает непрерывно, и каждая ячейка в качестве постоянной вероятности в единицу времени выполняет свою функцию обновления независимо от других ячеек.

$X_{ij}$$t_{ij}$$t_{ij}(0) \sim \operatorname{Exp}(\lambda)$$\lambda$ это параметр, значение которого можно выбрать произвольно.)

На каждом логическом временном шаге ячейки SCA обновляются следующим образом:

• $k, l$$i,j$$t_{kl}<t_{ij}$

• $X_{ij}$$X_{kl}$$\Delta t \sim \operatorname{Exp}(\lambda)$$t_{ij}$$t_{ij}+\Delta t$

Я полагаю, что это гарантирует, что ячейки будут обновлены в порядке, который может быть «декодирован», чтобы соответствовать исходному ACA, избегая коллизий и позволяя некоторым ячейкам обновляться параллельно. Однако из-за первого пункта выше, это означает, что большинство процессоров GPU будут в основном простаивать на каждом временном шаге SCA, что не идеально.

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

Я бы использовал первый вариант и использовал бы синхронный прогон AC перед (с использованием GPU), чтобы обнаружить коллизии, выполнить шаг шестиугольного AC, правилом которого является значение центральной ячейки = Sum (соседей), этот CA должен иметь семь состояний должны быть инициированы случайно выбранной ячейкой, и их состояние должно быть проверено перед запуском правила обновления для каждого графического процессора.

Пример 1. Значение соседней ячейки является общим

0 0 0 0 0 0 0 0

  0 0 1 0 0 0

0 0 0 0 0 0 0 0

  0 0 0 1 0 0

0 0 0 0 0 0 0 0

шаг CA, правилом которого является гексагональная центральная ячейка = Sum (соседи)

0 0 1 1 0 0 0

  0 1 1 1 0 0

0 0 1 2 1 0 0

  0 0 1 1 1 0

0 0 0 1 1 0 0

Пример 2. Значение ячейки для обновления учитывается как сосед с другой

0 0 0 0 0 0 0 0

  0 0 1 0 0 0

0 0 0 1 0 0 0

  0 0 0 0 0 0

0 0 0 0 0 0 0 0

После итерации

0 0 1 1 0 0 0

  0 1 2 2 0 0

0 0 2 2 1 0 0

  0 0 1 1 0 0

0 0 0 0 0 0 0 0

Образец 3. Нет связи

  0 0 0 0 0 0

0 0 1 0 0 0 0

  0 0 0 0 0 0

0 0 0 0 0 0 0 0

  0 0 0 1 0 0

0 0 0 0 0 0 0 0

После итерации

  0 1 1 0 0 0

0 1 1 1 0 0 0

  0 1 1 0 0 0

0 0 0 1 1 0 0

  0 0 1 1 1 0

0 0 0 1 1 0 0

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

Вы можете сделать это, используя атомарные операции, предусмотренные в CUDA, и массив, intсодержащий блокировки для каждой ячейки, например lock. Затем каждый поток выполняет следующее:

ci, cj = choose a pair at random.

int locked = 0;

/* Try to lock the cell ci. */
if ( atomicCAS( &lock[ci] , 0 , 1 ) == 0 ) {

    /* Try to lock the cell cj. */
    if ( atomicCAS( &lock[cj] , 0 , 1 ) == 0 ) {

        /* Now try to lock all the neigbourhood cells. */
        for ( cn = indices of all neighbours )
            if ( atomicCAS( &lock[cn] , 0 , 1 ) != 0 )
                break;

        /* If we hit a break above, we have to unroll all the locks. */
        if ( cn < number of neighbours ) {
            lock[ci] = 0;
            lock[cj] = 0;
            for ( int i = 0 ; i < cn ; i++ )
                lock[i] = 0;
            }

        /* Otherwise, we've successfully locked-down the neighbourhood. */
        else
            locked = 1;

        }

    /* Otherwise, back off. */
    else
        lock[ci] = 0;
    }

/* If we got everything locked-down... */
if ( locked ) {

    do whatever needs to be done...

    /* Release all the locks. */
    lock[ci] = 0;
    lock[cj] = 0;
    for ( int i = 0 ; i < cn ; i++ )
        lock[i] = 0;

    }

Обратите внимание, что этот подход, вероятно, не самый оптимальный, но он может обеспечить интересную отправную точку. Если между потоками существует много коллизий, то есть один или несколько на 32 нити (как в случае одной коллизии на деформацию), будет существенное отклонение ветвления. Кроме того, атомарные операции могут быть немного медленными, но поскольку вы выполняете только операции сравнения и замены, они должны масштабироваться нормально.

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

Заметьте также, что я быстр и свободен с обозначениями в петлях iнад соседями.

Приложение: Я был достаточно кавалерным, чтобы предположить, что вы можете просто отступить, когда пары сталкиваются. Если это не так, то вы можете заключить все, что и во второй строке в while-loop, и добавить breakв конце заключительного ifвыражения.

Все потоки должны будут ждать, пока последний не будет сделан, но если столкновения редки, вы должны быть в состоянии сойти с рук.

Приложение 2: Do не соблазниться добавить вызовы __syncthreads()в коде, особенно это зацикленной версия описана в предыдущем дополнении! Асинхронность имеет важное значение для предотвращения повторных столкновений в последнем случае.

Я ведущий разработчик LibGeoDecomp. Хотя я согласен с vanCompute, что вы можете эмулировать свой ACA с помощью CA, вы правы, что это не будет очень эффективным, так как только несколько ячеек на любом данном этапе должны быть обновлены. Это действительно очень интересное приложение - и с ним весело повозиться!

Я бы посоветовал вам объединить решения, предложенные jlopez1967 и Pedro: алгоритм Pedro хорошо улавливает параллелизм, но эти атомарные блокировки ужасно медленны. Решение jlopez1967 элегантно, когда дело доходит до обнаружения коллизий, но проверки всех nячеек, когда активен только меньший поднабор (теперь я буду предполагать, что есть некоторый параметр, kкоторый обозначает количество ячеек, которые должны обновляться одновременно), это явно непомерно.

__global__ void markPoints(Cell *grid, int gridWidth, int *posX, int *posY)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x, y;
    generateRandomCoord(&x, &y);
    posX[id] = x;
    posY[id] = y;
    grid[y * gridWidth + x].flag = 1;
}

__global__ void checkPoints(Cell *grid, int gridWidth, int *posX, int *posY, bool *active)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x = posX[id];
    int y = posY[id];
    int markedNeighbors = 
        grid[(y - 1) * gridWidth + x + 0].flag +
        grid[(y - 1) * gridWidth + x + 1].flag +
        grid[(y + 0) * gridWidth + x - 1].flag +
        grid[(y + 0) * gridWidth + x + 1].flag +
        grid[(y + 1) * gridWidth + x + 0].flag +
        grid[(y + 1) * gridWidth + x + 1].flag;
    active[id] = (markedNeighbors > 0);
}


__global__ void update(Cell *grid, int gridWidth, int *posX, int *posY, bool *active)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x = posX[id];
    int y = posY[id];
    grid[y * gridWidth + x].flag = 0;
    if (active[id]) {
        // do your fancy stuff here
    }
}

int main() 
{
  // alloc grid here, update up to k cells simultaneously
  int n = 1024 * 1024;
  int k = 1234;
  for (;;) {
      markPoints<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY);
      checkPoints<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY, active);
      update<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY, active);
  }
}

В отсутствие хорошей глобальной синхронизации на GPU, вам нужно вызвать несколько ядер для разных фаз. На Kepler от Nvidia вы можете переместить даже основной цикл на GPU, но я не ожидаю, что это сильно вырастет.

Алгоритмы достигают (настраиваемой) степени параллелизма. Я думаю, интересный вопрос - повлияют ли столкновения на ваше случайное распределение при увеличении k.

Я предлагаю вам увидеть эту ссылку http://www.wolfram.com/training/courses/hpc021.html примерно в 14:15 минут в видео, конечно, тренинг по математике, где они реализуют клеточные автоматы с использованием CUDA оттуда, и вы можете изменить его.