HPC - La mappa del gatto di Arnold

Moreno Marzolla moreno.marzolla@unibo.it

Ultimo aggiornamento: 2021-05-15

Anche la mappa del gatto di Arnold è una vecchia conoscenza, che abbiamo incontrato in una precedente esercitazione. In questo esercizio si chiede di realizzare un programma CUDA che trasforma una immagine mediante la mappa del gatto. Riportiamo nel seguito la descrizione del problema.

La mappa del gatto di Arnold è una funzione che trasforma una immagine \(P\) di dimensione \(N \times N\) in una nuova immagine \(P'\) delle stesse dimensioni. Per ogni \(0 \leq x < N,\ 0 \leq y < N\), il pixel di coordinate \((x,y)\) in \(P\) viene collocato nella posizione \((x',y')\) di \(P'\) dove:

\[ x' = (2x + y) \bmod N, \qquad y' = (x + y) \bmod N \]

("mod" è l'operatore modulo, corrispondente all'operatore % del linguaggio C). Si può assumere che le coordinate \((0, 0)\) indichino il pixel in alto a sinistra e le coordinate \((N-1, N-1)\) quello in basso a destra, in modo da poter indicizzare l'immagine come se fosse una matrice in linguaggio C. La Figura 1 mostra graficamente la trasformazione.

Figura 1: La mappa del gatto di Arnold
Figura 1: La mappa del gatto di Arnold

La mappa del gatto ha proprietà sorprendenti. Applicata ad una immagine ne produce una versione molto distorta. Applicata nuovamente a quest'ultima immagine, ne produce una ancora più distorta, e così via. Tuttavia, dopo un certo numero di iterazioni (il cui valore dipende da \(N\), ma che in ogni caso è sempre minore o uguale a \(3N\)) ricompare l'immagine di partenza! (Figura 2).

Figura 2: Alcune immagini ottenute iterando la mappa del gatto k volte
Figura 2: Alcune immagini ottenute iterando la mappa del gatto \(k\) volte

Il tempo minimo di ricorrenza per l'immagine cat1368.pgm di dimensione \(1368 \times 1368\) fornita come esempio è \(36\): iterando \(k\) volte della mappa del gatto si otterrà l'immagine originale se e solo se \(k\) è multiplo di 36. Non è nota alcuna formula analitica che leghi il tempo minimo di ricorrenza alla dimensione \(N\) dell'immagine.

Viene fornito un programma sequenziale che calcola la \(k\)-esima iterata della mappa del gatto usando la CPU. Il programma viene invocato specificando sulla riga di comando il numero di iterazioni \(k\). Il programma legge una immagine in formato PGM da standard input, e produce una nuova immagine su standard output ottenuta applicando \(k\) volte la mappa del gatto. Occorre ricordarsi di redirezionare lo standard output su un file, come indicato nelle istruzioni nel sorgente.

Per sfruttare il parallelismo offerto da CUDA è utile usare una griglia bidimensionale di thread block a loro volta bidimensionali, ciascuno con \(\mathit{BLKDIM} \times \mathit{BLKDIM}\) thread. Pertanto, data una immagine di \(N \times N\) pixel, sono necessari:

\[ (N + \mathit{BLKDIM} – 1) / \mathit{BLKDIM} \times (N + \mathit{BLKDIM} – 1) / \mathit{BLKDIM} \]

blocchi di dimensione \(\mathit{BLKDIM} \times \mathit{BLKDIM}\) per ricoprire interamente l'immagine. Ogni thread si occupa di calcolare una singola iterazione della mappa del gatto, copiando un pixel dell'immagine corrente nella posizione appropriata della nuova immagine. La segnatura del kernel sarà:

__global__ void cat_map_iter( unsigned char *cur, unsigned char *next, int N )

(dove \(N\) è la larghezza o altezza dell'immagine, che deve essere quadrata). Utilizzando il proprio ID e quello del blocco in cui si trova, ogni thread determina le coordinate \((x, y)\) del pixel su cui operare, e calcola le coordinate \((x', y')\) del pixel dopo l'applicazione di una iterazione della mappa del gatto. Per calcolare la \(k\)-esima iterata sarà quindi necessario invocare il kernel \(k\) volte, scambiando dopo ogni iterazione le immagini corrente e successiva come fatto dal programma seriale.

La soluzione precedente consente di parallelizzare il programma fornito apportando minime modifiche. Si può anche definire un kernel che calcoli direttamente la k-esima iterata della mappa del gatto con una singola invocazione. La segnatura del nuovo kernel sarà

__global__ void cat_map_iter_k( unsigned char *cur, unsigned char *next, int N, int k )

Come nel caso precedente, ciascun thread determina le coordinate \((x,y)\) del pixel di sua competenza. Le coordinate del pixel dopo \(k\) iterazioni si possono ottenere applicando lo schema seguente:

const int x = ...;
const int y = ...;
int xcur = x, ycur = y, xnext, ynext;

if ( x < N && y < N ) {
    while (k--) {
        xnext = (2*xcur + ycur) % N;
        ynext = (xcur + ycur) % N;
        xcur = xnext;
        ycur = ynext;
    }
    /* copia il pixel (x, y) dell'immagine corrente
    in posizione (xnext, ynext) della nuova immagine */
}

In questo modo è sufficiente una singola invocazione del kernel (anziché \(k\) come nel caso precedente) per ottenere l'immagine finale. Consiglio di misurare i tempi di esecuzione delle due alternative per capire se e di quanto la seconda soluzione è più efficiente della prima.

Per compilare:

    nvcc cuda-cat-map.cu -o cuda-cat-map

Per eseguire:

    ./cuda-cat-map k < input_file > output_file

Esempio:

    ./cuda-cat-map 100 < cat1368.pgm > cat1368.100.pgm

File