HPC - L'automa cellulare ANNEAL

Moreno Marzolla moreno.marzolla@unibo.it

Ultimo aggiornamento: 2021-11-28

In questo esercizio consideriamo un semplice automa cellulare binario in due dimensioni, denominato ANNEAL (noto anche come twisted majority rule). L'automa opera su un dominio quadrato di dimensione \(N \times N\), in cui ogni cella può avere valore 0 oppure 1. Si assume un dominio toroidale, in modo che ogni cella, incluse quelle sul bordo, abbia sempre otto celle adiacenti. Due celle si considerano adiacenti se hanno un lato oppure uno spigolo in comune.

L'automa evolve a istanti di tempo discreti \(t = 0, 1, 2, \ldots\). Lo stato di una cella al tempo \(t+1\) dipende dal proprio stato e da quello degli otto vicini al tempo \(t\). In dettaglio, per ogni cella \(x\) sia \(B_x\) il numero di celle con valore 1 presenti nell'intorno di dimensione \(3 \times 3\) centrato su \(x\) (si conta anche lo stato di \(x\), quindi si avrà sempre \(0 \leq B_x \leq 9\)). Se \(B_x=4\) oppure \(B_x \geq 6\), allora il nuovo stato della cella \(x\) è 1; in caso contrario il nuovo stato è 0. La Figura 1 mostra alcuni esempi.

Figura 1: Esempi di calcolo del nuovo stato della cella centrale di un blocco di dimensione 3 \times 3
Figura 1: Esempi di calcolo del nuovo stato della cella centrale di un blocco di dimensione \(3 \times 3\)

Come sempre in questi casi si devono utilizzare due griglie (domini) per rappresentare lo stato corrente dell'automa e lo stato al passo successivo. Lo stato delle celle viene sempre letto dalla griglia corrente, e i nuovi valori vengono sempre scritti nella griglia successiva. Quando il nuovo stato di tutte le celle è stato calcolato, si scambiano le griglie e si ripete.

Il dominio viene inizializzato ponendo ogni cella a 0 o 1 con uguale probabilità; di conseguenza, circa metà delle celle saranno nello stato 0 e l'altra metà sarà nello stato 1. La Figura 2 mostra l'evoluzione di una griglia di dimensione \(256 \times 256\) dopo 10, 100 e 1024 iterazioni. Si può osservare come le celle 0 e 1 tendano progressivamente ad addensarsi, pur con la presenza di piccole "bolle". È disponibile un breve video che mostra l'evoluzione dell'automa nel tempo.

Figura 2: Evoluzione dell'automa ANNEAL (video)
Figura 2: Evoluzione dell'automa ANNEAL (video)

Il file cuda-anneal.cu contiene una implementazione seriale dell'algoritmo che calcola e salva su un file l'evoluzione dopo \(K\) iterazioni dell'automa cellulare basato sulla regola ANNEAL. Scopo di questo esercizio è di modificare il programma per delegare alla GPU sia il calcolo del nuovo stato, sia la copia dei bordi del dominio (necessaria per simulare un dominio toroidale).

Alcuni suggerimenti:

Uso della shared memory

Questo programma potrebbe beneficiare dall'uso della memoria __shared__, dato che ogni cella del dominio viene letta 9 volte da 9 thread diversi. Tuttavia, sul server non si noterà alcun miglioramento delle prestazioni perché le GPU sono dotate di memoria cache e il numero di riletture non è abbastanza elevato da ammortizzare il costo della copia verso la memoria condivisa. Nonostante questo, è un esercizio istruttivo realizzare una versione che sfrutti la memoria shared.

Assumiamo che \(N\) sia un multiplo esatto di BLKDIM. Ciascun blocco di thread copia gli elementi della porzione di dominio di sua competenza in un buffer locale buf[BLKDIM+2][BLKDIM+2] che include due righe e due colonne (le prime e le ultime) di ghost cell, e calcolare il nuovo stato delle celle usando i dati nel buffer locale anziché accedendo alla memoria globale.

In situazioni del genere è utile usare due coppie di indici \((gi, gj)\) per indicare le posizioni delle celle nella matrice globale e \((li, lj)\) per le posizioni delle celle nel buffer locale. L'idea è che la cella di coordinate \((gi, gj)\) nella matrice globale corrisponda a quella di coordinate \((li, lj)\) nel buffer locale. Usando ghost cell sia a livello globale che a livello locale il calcolo delle coordinate può essere effettuato come segue:

    const int gi = 1 + threadIdx.y + blockIdx.y * blockDim.y;
    const int gj = 1 + threadIdx.x + blockIdx.x * blockDim.x;
    const int li = 1 + threadIdx.y;
    const int lj = 1 + threadIdx.x;
Figura 3: Copia dei dati dal dominio globale verso la shared memory
Figura 3: Copia dei dati dal dominio globale verso la shared memory

La parte più laboriosa è la copia dei dati dalla griglia globale al buffer locale. Usando \(\mathit{BLKDIM} \times \mathit{BLKDIM}\) thread per blocco, la copia della parte centrale (cioè tutto ad esclusione dell'area tratteggiata della Figura 3) si effettua con:

    buf[li][lj] = *IDX(cur, ext_n, gi, gj);

dove ext_n = (N + 2) è il lato del dominio inclusa la ghost area.

Figura 4: Thread attivi durante il riempimento del dominio locale
Figura 4: Thread attivi durante il riempimento del dominio locale

Per inizializzare la ghost area occorre procedere in tre fasi (Figura 4):

  1. La ghost area superiore e inferiore viene delegata ai thread della prima riga (quelli con \(li = 1\));

  2. La ghost area a sinistra e a destra viene delegata ai thread della prima colonna (quelli con \(lj = 1\));

  3. La ghost area negli angoli viene delegata al singolo thread con \((li, lj) = (1, 1)\).

(Si potrebbe essere tentati di collassare le fasi 1 e 2 in un'unica fase da far svolgere ad esempio ai thread della prima riga; questo sarebbe corretto, ma presenterebbe problemi nel caso si decidesse di generalizzare il codice a lati del dominio non necessariamente multipli di \(\mathit{BLKDIM}\)).

Si avrà in pratica la struttura seguente:

    if ( li == 1 ) {
        "riempi la cella buf[0][lj] e buf[BLKDIM+1][lj]"
    }
    if ( lj == 1 ) {
        "riempi la cella buf[li][0] e buf[li][BLKDIM+1]"
    }
    if ( li == 1 && lj == 1 ) {
        "riempi buf[0][0]"
        "riempi buf[0][BLKDIM+1]"
        "riempi buf[BLKDIM+1][0]"
        "riempi buf[BLKDIM+1][BLKDIM+1]"
    }

Chi vuole cimentarsi con una versione ancora più laboriosa può provare a modificare il codice per gestire anche il caso in cui la dimensione del dominio non sia multipla di BLKDIM. Prestare attenzione che non è sufficiente disattivare i thread al di fuori del dominio, ma bisogna modificare l'operazione di la copia della ghost area.

Per compilare senza usare shared memory:

    nvcc cuda-anneal.cu -o cuda-anneal

Per generare una immagine ad ogni passo:

    nvcc -DDUMPALL cuda-anneal.cu -o cuda-anneal

È possibile montare le immagini in un video in formato AVI/MPEG-4 con:

    ffmpeg -y -i "cuda-anneal-%05d.pbm" -vcodec mpeg4 cuda-anneal.avi

(il comando ffmpeg è già installato sul server).

Per compilare la soluzione abilitando la shared memory:

    nvcc -DUSE_SHARED cuda-anneal.cu -o cuda-anneal-shared

Per eseguire:

    ./cuda-anneal [steps [N]]

Esempio:

    ./cuda-anneal 64

File