Algorithmes parallèles (GPU) pour les automates cellulaires asynchrones

J'ai une collection de modèles informatiques qui pourraient être décrits comme des automates cellulaires asynchrones. Ces modèles ressemblent au modèle Ising, mais sont légèrement plus compliqués. Il semble que de tels modèles gagneraient à être exécutés sur un GPU plutôt que sur un CPU. Malheureusement, il n'est pas très simple de paralléliser un tel modèle, et je ne sais pas du tout comment procéder. Je suis conscient qu'il existe de la littérature sur le sujet, mais tout semble viser des informaticiens inconditionnels qui s'intéressent aux détails de la complexité algorithmique, plutôt que quelqu'un comme moi qui veut juste une description de quelque chose que je peux implémenter, et par conséquent, je le trouve plutôt intangible.

Pour plus de clarté, je ne recherche pas un algorithme optimal autant que quelque chose que je peux rapidement implémenter dans CUDA qui est susceptible de donner une accélération significative sur mon implémentation CPU. Le temps programmeur est beaucoup plus un facteur limitant que le temps ordinateur dans ce projet.

Je dois également préciser qu'un automate cellulaire asynchrone est quelque chose d'assez différent d'un automate synchrone, et les techniques de parallélisation des AC synchrones (comme la vie de Conway) ne peuvent pas être facilement adaptées à ce problème. La différence est qu'une autorité de certification synchrone met à jour chaque cellule simultanément à chaque pas de temps, tandis qu'une asynchrone met à jour une région locale choisie au hasard à chaque pas de temps, comme indiqué ci-dessous.

Les modèles que je souhaite paralléliser sont implémentés sur un réseau (généralement hexagonal) composé de ~ 100000 cellules (bien que j'aimerais en utiliser plus), et l'algorithme non parallélisé pour les exécuter ressemble à ceci:

1. Choisissez une paire de cellules voisine au hasard

2. Calculer une fonction "énergie" basée sur un voisinage local entourant ces cellules$\Delta E$

3. Avec une probabilité qui dépend de (avec β un paramètre), permutez les états des deux cellules ou ne faites rien.$e^{-\beta \Delta E}$$\beta$

4. Répétez les étapes ci-dessus indéfiniment.

Il y a aussi quelques complications liées aux conditions aux limites, mais j'imagine qu'elles ne poseront pas beaucoup de difficultés pour la parallélisation.

Il convient de mentionner que je m'intéresse à la dynamique transitoire de ces systèmes plutôt qu'à simplement l'état d'équilibre, j'ai donc besoin de quelque chose qui a une dynamique équivalente à ce qui précède, plutôt que simplement quelque chose qui s'approchera de la même distribution d'équilibre. (Les variations de l'algorithme du damier ne sont donc pas ce que je recherche.)

La principale difficulté pour paralléliser l'algorithme ci-dessus est les collisions. Étant donné que tous les calculs dépendent uniquement d'une région locale du réseau, il est possible que de nombreux sites de réseau soient mis à jour en parallèle, tant que leurs quartiers ne se chevauchent pas. La question est de savoir comment éviter de tels chevauchements. Je peux penser à plusieurs façons, mais je ne sais pas laquelle, s'il en est, est la meilleure à mettre en œuvre. Ce sont les suivants:

• Utilisez le processeur pour générer une liste de sites de grille aléatoires et vérifier les collisions. Lorsque le nombre de sites de grille est égal au nombre de processeurs GPU, ou si une collision est détectée, envoyez chaque ensemble de coordonnées à une unité GPU pour mettre à jour le site de grille correspondant. Cela serait facile à implémenter mais ne donnerait probablement pas beaucoup d'accélération, car la vérification des collisions sur le CPU ne serait probablement pas beaucoup moins chère que de faire toute la mise à jour sur le CPU.

• Divisez le réseau en régions (une par unité GPU) et demandez à une unité GPU de sélectionner et de mettre à jour au hasard les cellules de la grille dans sa région. Mais il y a beaucoup de problèmes avec cette idée que je ne sais pas comment résoudre, le plus évident étant ce qui devrait arriver exactement quand une unité choisit un quartier chevauchant le bord de sa région.

• Approximer le système comme suit: laisser le temps se dérouler par étapes discrètes. Divisez le réseau en un autreensemble de régions à chaque pas de temps selon un schéma prédéfini, et que chaque unité GPU sélectionne et met à jour au hasard une paire de cellules de grille dont le voisinage ne chevauche pas la frontière de la région. Comme les limites changent à chaque pas de temps, cette contrainte peut ne pas trop affecter la dynamique, tant que les régions sont relativement grandes. Cela semble facile à mettre en œuvre et susceptible d'être rapide, mais je ne sais pas dans quelle mesure il rapprochera la dynamique ou quel est le meilleur schéma pour choisir les limites de la région à chaque pas de temps. J'ai trouvé quelques références aux "automates cellulaires synchrones par blocs", qui peuvent ou non être les mêmes que cette idée. (Je ne sais pas car il semble que toutes les descriptions de la méthode soient en russe ou sont dans des sources auxquelles je n'ai pas accès.)

Mes questions spécifiques sont les suivantes:

• L'un des algorithmes ci-dessus est-il un moyen judicieux d'approcher la parallélisation GPU d'un modèle CA asynchrone?

• Y a-t-il une meilleure façon?

• Existe-t-il un code de bibliothèque pour ce type de problème?

• Où puis-je trouver une description claire en anglais de la méthode "block-synchronous"?

Le progrès

Je pense avoir trouvé un moyen de paralléliser une autorité de certification asynchrone qui pourrait convenir. L'algorithme décrit ci-dessous est pour une autorité de certification asynchrone normale qui met à jour une seule cellule à la fois, plutôt qu'une paire de cellules voisine comme la mienne. Il y a quelques problèmes à généraliser à mon cas spécifique, mais je pense avoir une idée de la façon de les résoudre. Cependant, je ne sais pas quel avantage de vitesse cela donnera, pour les raisons décrites ci-dessous.

L'idée est de remplacer l'AC asynchrone (désormais ACA) par une AC synchrone stochastique (SCA) qui se comporte de manière équivalente. Pour ce faire, nous imaginons d'abord que l'ACA est un processus de Poisson. Autrement dit, le temps se déroule en continu, et chaque cellule comme une probabilité constante par unité de temps d'exécuter sa fonction de mise à jour, indépendamment des autres cellules.

$X_{ij}$$t_{ij}$$t_{ij}(0) \sim \operatorname{Exp}(\lambda)$$\lambda$ est un paramètre dont la valeur peut être choisie arbitrairement.)

À chaque pas de temps logique, les cellules du SCA sont mises à jour comme suit:

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

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

Je crois que cela garantit que les cellules seront mises à jour dans un ordre qui peut être "décodé" pour correspondre à l'ACA d'origine, tout en évitant les collisions et en permettant à certaines cellules d'être mises à jour en parallèle. Cependant, en raison du premier point ci-dessus, cela signifie que la plupart des processeurs GPU seront principalement inactifs à chaque pas de temps du SCA, ce qui est loin d'être idéal.

Je dois réfléchir un peu plus à la question de savoir si les performances de cet algorithme peuvent être améliorées et comment étendre cet algorithme pour faire face au cas où plusieurs cellules sont mises à jour simultanément dans l'ACA. Cependant, cela semble prometteur, j'ai donc pensé le décrire ici au cas où quelqu'un (a) connaîtrait quelque chose de similaire dans la littérature, ou (b) pourrait offrir un aperçu de ces problèmes restants.

J'utiliserais la première option et utiliserais un AC synchrone avant (en utilisant le GPU), pour détecter les collisions, exécuter une étape d'un AC hexagonal dont la règle est la valeur de la cellule centrale = Sum (voisins), cette autorité de certification doit avoir sept états doivent être initiés avec une cellule sélectionnée au hasard, et leur état vérifié avant d'exécuter la règle de mise à jour pour chaque processeur.

Exemple 1. La valeur d'une cellule voisine est partagée

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

une étape d'une autorité de certification dont la règle est la cellule centrale hexagonale = somme (voisins)

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

Exemple 2. La valeur d'une cellule à mettre à jour est prise en compte comme voisine de l'autre

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

Après l'itération

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

Échantillon 3. Il n'y a pas de relation

  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

Après l'itération

  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

Suite à vos réponses à mes questions dans les commentaires ci-dessus, je vous suggère d'essayer une approche basée sur les verrous dans laquelle chaque thread essaie de verrouiller le quartier qu'il mettra à jour avant de calculer la mise à jour réelle.

Vous pouvez le faire en utilisant les opérations atomiques prévues dans CUDA et un tableau intcontenant les verrous pour chaque cellule, par exemple lock. Chaque thread effectue ensuite les opérations suivantes:

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;

    }

Notez que cette approche n'est probablement pas la plus optimale, mais elle pourrait fournir un point de départ intéressant. S'il y a beaucoup de collisions entre les threads, c'est-à-dire une ou plusieurs pour 32 threads (comme dans une collision par chaîne), il y aura alors un bon détournement de branche. De plus, les opérations atomiques peuvent être un peu lentes, mais comme vous ne faites que des opérations de comparaison et d'échange, elles devraient évoluer correctement.

La surcharge de verrouillage peut sembler intimidante, mais ce n'est vraiment que quelques affectations et branches, pas beaucoup plus.

Notez également que je suis rapide et lâche avec la notation dans les boucles de isur les voisins.

Addendum: j'étais assez cavalière pour supposer que vous pouviez simplement reculer quand les paires entrent en collision. Si ce n'est pas le cas, vous pouvez tout envelopper à partir de la deuxième ligne dans une whileboucle et ajouter un breakà la fin de la ifdéclaration finale .

Tous les threads devront alors attendre que le dernier soit terminé, mais si les collisions sont rares, vous devriez pouvoir vous en tirer.

Addendum 2: Ne soyez pas tenté d'ajouter des appels à __syncthreads()n'importe où dans ce code, en particulier à la version en boucle décrite dans l'addendum précédent! L'asynchronicité est essentielle pour éviter les collisions répétées dans ce dernier cas.

Je suis le développeur principal de LibGeoDecomp. Bien que je convienne avec vanCompute que vous pouvez émuler votre ACA avec une autorité de certification, vous avez raison de dire que cela ne serait pas très efficace, car seules quelques cellules à une étape donnée sont censées être mises à jour. Il s'agit en effet d'une application très intéressante - et amusante à bricoler!

Je vous suggère de combiner les solutions proposées par jlopez1967 et Pedro: l'algorithme de Pedro capture bien le parallélisme, mais ces verrous atomiques sont terriblement lents. La solution de jlopez1967 est élégante en ce qui concerne la détection des collisions, mais la vérification de toutes les ncellules, quand seul un sous-ensemble plus petit (je suppose désormais qu'il existe un paramètre kqui dénote le nombre de cellules à mettre à jour simultanément) est actif, est clairement prohibitif.

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

En l'absence d'une bonne synchronisation globale sur le GPU, vous devez invoquer plusieurs noyaux pour les différentes phases. Sur Kepler de Nvidia, vous pouvez déplacer même la boucle principale vers le GPU, mais je ne m'attends pas à ce que cela gagne beaucoup.

Les algorithmes atteignent un degré (configurable) de parallélisme. Je suppose que la question intéressante est de savoir si les collisions affecteront votre distribution aléatoire lorsque vous augmentez k.

Je vous suggère que vous voyez ce lien http://www.wolfram.com/training/courses/hpc021.html environ 14:15 minutes dans la vidéo bien sûr, une formation mathématique où ils font une implémentation d'un automate cellulaire en utilisant CUDA , à partir de là et vous pouvez le modifier.