Algorytmy równoległe (GPU) dla asynchronicznych automatów komórkowych

Mam kolekcję modeli obliczeniowych, które można opisać jako asynchroniczne automaty komórkowe. Modele te przypominają model Isinga, ale są nieco bardziej skomplikowane. Wydaje się, że takie modele skorzystałyby na GPU, a nie na CPU. Niestety równoległość takiego modelu nie jest łatwa i wcale nie jest dla mnie jasne, jak sobie z tym poradzić. Wiem, że istnieje literatura na ten temat, ale wszystko to wydaje się być skierowane do zapalonych informatyków, którzy są zainteresowani szczegółami złożoności algorytmicznej, a nie do kogoś takiego jak ja, który chce tylko opisu czegoś, co mogę zaimplementować, i w konsekwencji uważam, że jest to raczej nieprzeniknione.

Dla jasności nie szukam optymalnego algorytmu, a nie czegoś, co mogę szybko zaimplementować w CUDA, co może znacznie przyspieszyć moją implementację procesora. W tym projekcie czas programisty jest znacznie bardziej ograniczającym czynnikiem niż czas komputerowy.

Powinienem również wyjaśnić, że asynchroniczny automat komórkowy jest czymś innym niż synchroniczny, a techniki równoległego synchronicznego CA (takie jak życie Conwaya) nie mogą być łatwo dostosowane do tego problemu. Różnica polega na tym, że synchroniczny urząd certyfikacji aktualizuje każdą komórkę jednocześnie na każdym etapie, podczas gdy asynchroniczny aktualizuje losowo wybrany region lokalny na każdym etapie, jak opisano poniżej.

Modele, które chcę zrównoleglać, są zaimplementowane w sieci (zwykle heksagonalnej) składającej się z ~ 100000 komórek (choć chciałbym użyć ich więcej), a niesparalizowany algorytm ich uruchamiania wygląda następująco:

1. Wybierz losowo sąsiednią parę komórek

2. Oblicz funkcję „energii” na podstawie lokalnego sąsiedztwa otaczającego te komórki$\Delta E$

3. Z prawdopodobieństwem zależnym od (z parametrem a) albo zamień stany dwóch komórek, albo nic nie rób.$e^{-\beta \Delta E}$$\beta$

4. Powtórz powyższe kroki w nieskończoność.

Istnieją również pewne komplikacje związane z warunkami brzegowymi, ale wyobrażam sobie, że nie będą one stanowić dużego problemu dla równoległości.

Warto wspomnieć, że interesuje mnie przejściowa dynamika tych układów, a nie tylko stan równowagi, więc potrzebuję czegoś, co ma równoważną dynamikę do powyższych, a nie tylko czegoś, co zbliży się do tego samego rozkładu równowagi. (Więc odmiany algorytmu szachownicy nie są tym, czego szukam.)

Główną trudnością w równoległości powyższego algorytmu są kolizje. Ponieważ wszystkie obliczenia zależą tylko od lokalnego regionu sieci, wiele witryn sieci może być aktualizowanych równolegle, o ile ich sąsiedztwo się nie nakłada. Pytanie brzmi, jak uniknąć takiego nakładania się. Mogę wymyślić kilka sposobów, ale nie wiem, który z nich najlepiej wdrożyć. Są to:

• Użyj procesora, aby wygenerować listę losowych witryn sieci i sprawdzić kolizje. Gdy liczba miejsc siatki jest równa liczbie procesorów GPU lub w przypadku wykrycia kolizji, wyślij każdy zestaw współrzędnych do jednostki GPU, aby zaktualizować odpowiednie miejsce sieci. Byłoby to łatwe do wdrożenia, ale prawdopodobnie nie dałoby większego przyspieszenia, ponieważ sprawdzanie kolizji na CPU prawdopodobnie nie byłoby o wiele tańsze niż wykonanie całej aktualizacji na CPU.

• Podziel sieć na regiony (po jednej na jednostkę GPU) i niech jedna jednostka GPU odpowiada za losowy wybór i aktualizację komórek siatki w swoim regionie. Istnieje jednak wiele problemów z tym pomysłem, których nie umiem rozwiązać, najbardziej oczywistym jest to, co dokładnie powinno się zdarzyć, gdy jednostka wybierze dzielnicę pokrywającą się z krawędzią swojego regionu.

• Przybliż system w następujący sposób: pozwól, aby czas przebiegał w dyskretnych krokach. Podziel sieć na innązestaw regionów na każdym etapie za każdym razem zgodnie z pewnym wstępnie zdefiniowanym schematem, a każda jednostka GPU losowo wybiera i aktualizuje parę komórek siatki, których sąsiedztwo nie zachodzi na granicę regionu. Ponieważ granice zmieniają się za każdym razem, ograniczenie to może nie wpływać zbytnio na dynamikę, o ile regiony są stosunkowo duże. Wydaje się to łatwe do wdrożenia i prawdopodobnie szybkie, ale nie wiem, jak dobrze będzie to przybliżać dynamikę, ani jaki jest najlepszy schemat wyboru granic regionu na każdym kroku. Znalazłem kilka odniesień do „blokowo-synchronicznych automatów komórkowych”, które mogą, ale nie muszą być takie same jak ten pomysł. (Nie wiem, ponieważ wydaje się, że wszystkie opisy metody są albo w języku rosyjskim, albo w źródłach, do których nie mam dostępu.)

Moje konkretne pytania są następujące:

• Czy którykolwiek z powyższych algorytmów jest rozsądnym sposobem podejścia do równoległości GPU asynchronicznego modelu CA?

• Czy jest lepszy sposób?

• Czy istnieje kod biblioteki dla tego typu problemu?

• Gdzie mogę znaleźć jasny anglojęzyczny opis metody „blokowo-synchronicznej”?

Postęp

Wydaje mi się, że wpadłem na pomysł równoległego asynchronicznego urzędu certyfikacji, który może być odpowiedni. Algorytm przedstawiony poniżej dotyczy normalnego asynchronicznego urzędu certyfikacji, który aktualizuje tylko jedną komórkę na raz, a nie sąsiednią parę komórek, tak jak moja. Istnieją pewne problemy z ogólnym uogólnieniem tego przypadku, ale myślę, że mam pomysł, jak je rozwiązać. Jednak nie jestem pewien, ile korzyści z prędkości przyniesie, z powodów omówionych poniżej.

Chodzi o to, aby zastąpić asynchroniczny CA (odtąd ACA) stochastycznym synchronicznym CA (SCA), który zachowuje się równorzędnie. Aby to zrobić, najpierw wyobrażamy sobie, że ACA jest procesem Poissona. Oznacza to, że czas płynie w sposób ciągły, a każda komórka jako stałe prawdopodobieństwo wykonania funkcji aktualizacji na jednostkę czasu, niezależnie od innych komórek.

Konstruujemy SCA, którego komórki przechowują dwie rzeczy: stan komórki (tj. Dane, które będą przechowywane w każdej komórce w implementacji sekwencyjnej) oraz liczbę zmiennoprzecinkową reprezentującą (ciągły ) godzina, o której nastąpi kolejna aktualizacja. Ten ciągły czas nie odpowiada etapom aktualizacji SCA. Odniosę się do tego drugiego jako do „czasu logicznego”. Wartości czasu są inicjowane losowo zgodnie z rozkładem wykładniczym: . (Gdzie jest parametrem, którego wartość można wybrać dowolnie.) t i j t i j ( 0 ) ∼ Exp ( λ ) $X_{ij}$$t_{ij}$$t_{ij}(0) \sim \operatorname{Exp}(\lambda)$$\lambda$

Na każdym logicznym etapie czasu komórki SCA są aktualizowane w następujący sposób:

• Jeśli dla dowolnego sąsiedztwie , czas , nie rób nic.i , j t k l < t i $k, l$$i,j$$t_{kl}<t_{ij}$

• W przeciwnym razie (1) zaktualizuj stan zgodnie ze stanami sąsiednich komórek, stosując tę ​​samą regułę, co oryginalny ACA; i (2) wygeneruj losową wartość nazwa i zaktualizuj do . X k l Δ t ∼ Exp ( λ ) t i j t i j + Δ $X_{ij}$$X_{kl}$$\Delta t \sim \operatorname{Exp}(\lambda)$$t_{ij}$$t_{ij}+\Delta t$

Wierzę, że to gwarantuje, że komórki zostaną zaktualizowane w kolejności, którą można „zdekodować”, aby odpowiadały pierwotnemu ACA, unikając kolizji i umożliwiając równoległą aktualizację niektórych komórek. Jednak ze względu na pierwszy punkt powyżej oznacza to, że większość procesorów GPU będzie w większości bezczynna na każdym kroku SCA, co jest mniej niż idealne.

Muszę się zastanowić, czy wydajność tego algorytmu można poprawić, i jak rozszerzyć ten algorytm, aby poradził sobie z przypadkiem, gdy wiele komórek jest aktualizowanych jednocześnie w ACA. Wygląda jednak obiecująco, więc pomyślałem, że opisałbym to tutaj na wypadek, gdyby ktokolwiek (a) wiedział o czymś podobnym w literaturze lub (b) mógł zaoferować wgląd w pozostałe kwestie.

Użyłbym pierwszej opcji i użyłbym synchronicznego przebiegu prądu przemiennego przed (za pomocą GPU), aby wykryć kolizje, wykonać krok sześciokątnego prądu przemiennego, którego regułą jest wartość komórki środkowej = Suma (sąsiedzi), ten urząd certyfikacji musi mieć siedem stanów powinno zostać zainicjowanych losowo wybraną komórką, a ich status zweryfikowany przed uruchomieniem reguły aktualizacji dla każdej GPU.

Próbka 1. Wartość sąsiedniej komórki jest wspólna

0 0 0 0 0 0 0

  0 0 1 0 0 0

0 0 0 0 0 0 0

  0 0 0 1 0 0

0 0 0 0 0 0 0

krok CA, którego regułą jest heksagonalna komórka centralna = Suma (sąsiedzi)

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

Próbka 2. Wartość komórki do aktualizacji jest brana pod uwagę jako sąsiad z drugiej

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

Po iteracji

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

Próbka 3. Brak związku

  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 1 0 0

0 0 0 0 0 0 0

Po iteracji

  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

Po udzieleniu odpowiedzi na moje pytania w powyższych komentarzach sugeruję wypróbowanie metody blokowania, w której każdy wątek próbuje zablokować okolicę, którą będzie aktualizować, przed obliczeniem rzeczywistej aktualizacji.

Możesz to zrobić za pomocą operacji atomowych przewidzianych w CUDA i tablicy intzawierającej blokady dla każdej komórki, np lock. Każdy wątek wykonuje następujące czynności:

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;

    }

Zauważ, że to podejście nie jest prawdopodobnie najbardziej optymalne, ale może stanowić ciekawy punkt wyjścia. Jeśli występuje wiele kolizji między wątkami, tj. Jeden lub więcej na 32 wątki (jak w jednym kolizji na osnowę), nastąpi spore przesunięcie gałęzi. Również operacje atomowe mogą być nieco powolne, ale ponieważ wykonujesz tylko operacje porównywania i zamiany, powinno się skalować.

Blokowanie może wydawać się zastraszające, ale tak naprawdę to tylko kilka zadań i gałęzi, niewiele więcej.

Zauważ też, że jestem szybki i swobodny z zapisem w pętlach inad sąsiadami.

Dodatek: Byłem na tyle kawalerski, aby założyć, że można po prostu wycofać się, kiedy pary się zderzą. Jeśli tak nie jest, możesz owinąć wszystko od drugiej linii w whilepętlę i dodać znak breakna końcu ostatniej ifwypowiedzi.

Wszystkie wątki będą musiały poczekać, aż ostatni zostanie zrobiony, ale jeśli kolizje są rzadkie, powinieneś być w stanie uciec.

Dodatek 2: Czy nie ulec pokusie, aby dodać do połączenia __syncthreads()gdziekolwiek w tym kodzie, zwłaszcza, że jest to wersja pętli opisany w poprzednim dodatku! Asynchroniczność jest niezbędna w celu uniknięcia powtarzających się kolizji w tym drugim przypadku.

Jestem głównym programistą LibGeoDecomp. Chociaż zgadzam się z vanCompute, że możesz emulować swoją ACA za pomocą urzędu certyfikacji, masz rację, że nie byłoby to bardzo wydajne, ponieważ tylko kilka komórek w danym kroku ma zostać zaktualizowanych. Jest to rzeczywiście bardzo interesująca aplikacja - i majstrowanie przy niej jest zabawne!

Sugeruję połączenie rozwiązań zaproponowanych przez jlopez1967 i Pedro: algorytm Pedro dobrze przechwytuje równoległość, ale te blokady atomowe są strasznie powolne. Rozwiązanie jlopez1967 jest eleganckie, jeśli chodzi o wykrywanie kolizji, ale sprawdzanie wszystkich nkomórek, gdy aktywny jest tylko mniejszy podzbiór (od tej pory założę, że istnieje parametr kokreślający liczbę komórek do aktualizacji jednocześnie), jest wyraźnie wygórowany.

__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);
  }
}

W przypadku braku dobrej globalnej synchronizacji na GPU musisz wywołać wiele jąder dla różnych faz. W Kepler Nvidii możesz przenieść nawet główną pętlę do GPU, ale nie spodziewam się, że to wiele zyska.

Algorytmy osiągają (konfigurowalny) stopień równoległości. Wydaje mi się, że interesujące pytanie dotyczy tego, czy kolizje wpłyną na losowy rozkład, gdy wzrośnie k.

Sugeruję, abyś zobaczył ten link http://www.wolfram.com/training/courses/hpc021.html około 14:15 minut do filmu, oczywiście, szkolenia matematyczne, w których dokonują implementacji automatów komórkowych przy użyciu CUDA , stamtąd i możesz go zmodyfikować.