HPC - La mappa del gatto di Arnold

Scopo di questo esercizio è sviluppare una versione SIMD di una funzione che calcola l’iterata della mappa del gatto di Arnold, una vecchia conoscenza che abbiamo già incontrato in altre esercitazioni. Riportiamo nel seguito la descrizione del problema.

La mappa del gatto 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:

(“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.

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! (si veda la Figura 2).

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. La struttura della funzione che calcola la k-esima iterata della mappa del gatto è molto semplice:

for (y=0; y<N; y++) {
    for (x=0; x<N; x++) {
        /* calcola le coordinate (xnew, ynew) del punto (x, y)
            dopo k applicazioni della mappa del gatto */
        next[xnew + ynew*N] = cur[x+y*N];
    }
}

Per sfruttare il parallelismo SIMD possiamo ragionare come segue: anziché calcolare le nuove coordinate di un punto alla volta, calcoliamo le coordinate di quattro punti adiacenti \((x, y)\), \((x+1,y)\), \((x+2,y)\), \((x+3,y)\) usando i vector datatype del compilatore. Per fare questo, definiamo le seguenti variabili di tipo v4i (vettori SIMD di 4 interi):

    typedef int v4i __attribute__((vector_size(16)));
    #define VLEN (sizeof(v4i)/sizeof(int))

Posto \(vx = \{x, x+1, x+2, x+3\}\), \(vy = \{y, y, y, y\}\), possiamo applicare ad essi le stesse operazioni aritmetiche applicate agli scalari x e y per ottenere le nuove coordinate vxnew, vynew. Fatto questo, al posto della singola istruzione:

    next[xnew + ynew*N] = cur[x+y*N];

per spostare materialmente i pixel nella nuova posizione occorre eseguire quattro istruzioni scalari:

    next[vxnew[0] + vynew[0]*N] = cur[vx[0] + vy[0]*N];
    next[vxnew[1] + vynew[1]*N] = cur[vx[1] + vy[1]*N];
    next[vxnew[2] + vynew[2]*N] = cur[vx[2] + vy[2]*N];
    next[vxnew[3] + vynew[3]*N] = cur[vx[3] + vy[3]*n];

Estensione

Le prestazioni della versione SIMD dell’algoritmo della mappa del gatto dovrebbero risultare solo marginalmente migliori della versione scalare (potrebbero addirittura essere peggiori). Analizzando il codice assembly prodotto dal compilatore, si scopre che il calcolo del modulo nelle due espressioni

    vxnew = (2*vxold+vyold) % N;
    vynew = (vxold + vyold) % N;

viene realizzato usando operazioni scalari. Consultando la lista dei SIMD intrinsics sul sito di Intel si scopre che non esiste una istruzione SIMD che realizzi la divisione intera. Per migliorare le prestazioni del programma occorre quindi ingegnarsi per calcolare i moduli senza fare uso della divisione. Ragionando in termini scalari, osserviamo che se \(0 \leq xold < N\) e \(0 \leq yold < N\), allora si ha necessariamente che \(0 \leq 2 \times xold + yold < 3N\) e \(0 \leq xold+yold < 2N\).

Pertanto, sempre in termini scalari, possiamo realizzare il calcolo di xnew e ynew come segue:

    xnew = (2*xold + yold);
    if (xnew >= N) { xnew = xnew - N; }
    if (xnew >= N) { xnew = xnew - N; }
    ynew = (xold + yold);
    if (ynew >= N) { ynew = ynew - N; }

Il codice precedente è meno leggibile della versione che usa l’operatore modulo, ma ha il vantaggio di poter essere vettorizzato ricorrendo al meccanismo di “selection and masking” visto a lezione. Ad esempio, l’istruzione

    if (xnew >= N) { xnew = xnew - N; }

    const v4i mask = (xnew >= N);
    xnew = (mask & (xnew - N)) | (mask & xnew);

    const v4i mask = (xnew >= N);
    xnew = xnew - (mask & N);

Si ottiene in questo modo un programma più complesso della versione scalare, ma più veloce in quanto si riesce a sfruttare al meglio le istruzioni SIMD.

Estensione

File