HPC - Automa cellulare della "regola 30"

Questo esercizio chiede di implementare mediante CUDA l'automa cellulare della regola 30 che abbiamo già incontrato in una precedente esercitazione. Per comodità riportiamo la descrizione del funzionamento dell'automa.

L'automa è costituito da un array \(x[N]\) di \(N\) interi, ciascuno dei quali può avere valore 0 oppure 1. Lo stato dell'automa evolve durante istanti discreti nel tempo: lo stato di una cella al tempo \(t\) dipende dal suo stato e da quello dei due vicini al tempo \(t-1\). Assumiamo un dominio ciclico, per cui i vicini della cella \(x[0]\) sono \(x[N-1]\) e \(x[1]\) e i vicini della cella \(x[N-1]\) sono \(x[N-2]\) e \(x[0]\).

Dati i valori correnti \(pqr\) di tre celle adiacenti, il nuovo valore \(q'\) della cella centrale è determinato in base alla tabella seguente (■ = 1, □ = 0):

(si noti che la sequenza □□□■■■■□ = 00011110, che si legge sulla seconda riga della tabella, rappresenta il numero 30 in binario, da cui il nome "regola 30").

Il file cuda-rule30.cu contiene lo scheletro di una implementazione seriale dell'algoritmo che calcola l'evoluzione dell'automa. Inizialmente tutte le celle sono nello stato 0, ad eccezione di quella in posizione \(N/2\) che è nello stato 1. Il programma accetta sulla riga di comando la dimensione del dominio \(N\) e il numero di passi da calcolare (se non indicati, si usano dei valori di default). Al termine dell'esecuzione il programma produce il file cuda-rule30.pbm contenente una immagine simile alla Figura 1.

Ogni riga dell'immagine rappresenta lo stato dell'automa in un dato istante di tempo; il colore di ogni pixel indica lo stato di ciascuna cella (1 = nero, 0 = bianco). Il tempo avanza dall'alto verso il basso, quindi la prima riga indica lo stato al tempo 0, la seconda riga lo stato al tempo 1 e così via.

Scopo di questo esercizio è di svilupparne una versione parallela in cui il calcolo dei nuovi stati (cioè di ogni riga dell'immagine) venga realizzato da thread CUDA. In particolare, la funzione rule30() dovrà essere trasformata in un kernel che viene invocato per calcolare un passo di evoluzione dell'automa. Assumere che \(N\) sia multiplo del numero di thread per blocco (BLKDIM).

Si consiglia di iniziare con una versione in cui si opera direttamente sulla memoria globale senza usare memoria condivisa (__shared__); questa prima versione si può ricavare molto velocemente partendo dal codice seriale fornito come esempio.

Da quanto detto a lezione, l'uso della memoria condivisa può essere utile perché ogni valore di stato nella memoria globale viene letto più volte (tre, per la precisione). Realizzare una seconda versione del programma che sfrutti la memoria condivisa. L'idea è la stessa dell'esempio visto a lezione relativo ad una computazione di tipo stencil, con la differenza che in questo caso il raggio della stencil vale 1, cioè il nuovo stato di ogni cella dipende dallo stato precedente di quella cella e dei due vicini. Si presti però attenzione che, a differenza dell'esempio visto a lezione, qui si assume un dominio ciclico. La parte più delicata sarà quindi la copia dei valori dalla memoria globale alla memoria condivisa.

Aiutandoci con la Figura 2 procediamo come segue:

• d_cur[] è la copia nella memoria del device dello stato corrente dell'automa.

• Definire un kernel per riempire le ghost cell di d_cur[]; tale kernel, chiamato ad esempio fill_ghost(...), deve essere eseguito da un singolo thread e quindi andrà invocato come fill_ghost<<<1, 1>>>(...)

• Definire un secondo kernel 1D per aggiornare il dominio. Ciascun thread block definisce un array __shared__ chiamato ad esempio buf[BLKDIM+2]; usiamo BLKDIM+2 elementi perché dobbiamo includere una ghost cell a sinistra e una a destra per calcolare lo stato successivo dei BLKDIM elementi centrali senza dover ricorrere a ulteriori letture dalla memoria globale;

• Ciascun thread determina l'indice lindex dell'elemento locale (nell'array buf[]), e l'indice gindex dell'elemento globale (nell'array cur[] in memoria globale) su cui deve operare. Estendiamo l'array cur[] con due ghost cell (una per estremità), come pure l'array buf[] in memoria condivisa. Quindi gli indici vanno calcolati come:

    const int lindex = 1 + threadIdx.x;
    const int gindex = 1 + threadIdx.x + blockIdx.x * blockDim.x;

• Ciascun thread copia un elemento dalla memoria globale alla memoria condivisa:

    buf[lindex] = cur[gindex];

• Il primo thread di ciascun blocco inizializza le ghost cell di buf[] ad esempio con:

    if (0 == threadIdx.x) {
        buf[0] = cur[gindex-1];
        buf[BLKDIM + 1] = cur[gindex + BLKDIM];
    }

Per generare l'evoluzione dell'automa, le nuove configurazioni vanno trasferite dal device all'host al termine di ogni esecuzione del kernel.

Per compilare:

    nvcc cuda-rule30.cu -o cuda-rule30

Per eseguire:

    ./cuda-rule30 [width [steps]]

Esempio:

    ./cuda-rule30 1024 1024

Il risultato sarà nel file cuda-rule30.pbm

File

• cuda-rule30.cu
• hpc.h