GPGPU.IT

Nella prima parte di questo corso abbiamo visto una prima generale descrizione del modello architetturale di CUDA e realizzato il primo programmino in grado di allocare strutture dati nella memoria della GPU e spostare i dati tra di esse. In questa seconda puntata comprenderemo cosa significa definire una configurazione di esecuzione in termini di blocchi e thread e realizzeremo il primo kernel in grado di eseguire effettivamente dei calcoli.

Qualificatori di funzione

Un kernel CUDA è una particolare funzione C che, invocata dall'host (CPU), viene eseguita sul device (la GPU). I kernel hanno delle caratteristiche ben precise. Innazitutto, devono avere void come tipo di ritorno (cioè, non devono ritornare valori), non possono essere ricorsive, non possono avere un numero di parametri variabile e non possono usare variabili di tipo statico. Inoltre, in generale, è bene considerare che possano accedere solo alla memoria della GPU, dunque a strutture dati che siano allocate nella global memory del device. A partire dalla versione 2.2 di CUDA questo non è sempre vero, visto che è possibile sfruttare il meccanismo di Zero Copy con la primitiva cudaMallocHost per consentire alla GPU di leggere e scrivere dati direttamente nella memoria dell'Host, ma vedremo solo successivamente quali siano le sue implicazioni.

Nel C for CUDA le funzioni facenti parte di un kernel vanno identificate tramite appositi qualificatori da anteporre alla loro dichiarazione. Esistono tre qualificatori:

• __global__ : identifica la funzione invocata dal codice dell'host, è effettivamente la funzione principale del kernel per la quale valgono le caratteristiche sopra elencate
• __device__: identifica una funzione chiamata dal codice che gira sulla GPU, si comporta come una normale funzione C (dunque non valgono le limitazioni per i kernel), ma non può essere invocata dall'host
• __host__: è una funzione chiamata dall'host, si comporta come una normale funzione C,  ma viene eseguita solo sull'host.

Le funzioni __device__ ed __host__ possono essere combinate, per esempio per generare codice che funzioni sia sulla GPU che sulla CPU, sfruttando il meccanismo di overloading degli operatori.

L'invocazione di un kernel avviene all'interno del codice eseguito sulla CPU, chiamandone il nome e specificando la configurazione di esecuzione all'interno di appositi identificatori (<<< e >>>) con la seguente sintassi:

Kernel<<<dim3 grid, dim3 block>>>(parameter1, paramenter2, ...);

grid e block sono delle variabili che specificano il numero di blocchi in cui è diviso il kernel e il numero di thread all'interno di ciascun blocco. Volendo, possono essere direttamente sostituite con dei numeri.

Configurazione di esecuzione

dim3 è uno dei tipi predefiniti di CUDA. Come il nome lascia intuire, è un tipo di vettore tridimensionale che permette di specificare dunque fino a tre valori per la variabile dichiarata. Valori eventualmente non dichiarati vengono automaticamente posti a 1.

Come abbiamo precedentemente detto, un kernel dal punto di vista del codice è anche definito grid (griglia). La griglia è composta da un insieme di blocchi, che può essere bidimensionale (2D). I blocchi sono infine costituiti dai thread, che possono essere organizzati tridimensionalmente. Dunque, specificare una esecuzione di configurazione, significa dichiarare ad esempio le seguenti variabili:

• dim3 grid(16, 16);
• dim3 block(16,16);

In questa maniera, impostiamo una variabile grid e una variabile block inizializzate a 16, 16, 1. Si può così lanciare il kernel:

Kernel<<<grid, block>>>( parameter1);

è bene ricordare che un kernel può essere lanciato anche senza bisogno di usare le variabili di tipo dim3 e specificando direttamente i valori all'interno dell'invocazione:

Kernel<<<16,256>>>(parameter1);

Questo comando permette di lanciare un kernel composto da 16 blocchi con 256 thread al suo interno. Gli altri valori vengono automaticamente posti a 1. In generale, però, consigliamo di usare codice il più parametrizzato possibile, in modo da favorire la sperimentazione di differenti configurazioni di esecuzione per trovare la migliore possibile. Vedremo in seguito che c'è la possibilità di cercare di individuare una buona configurazione di esecuzione, almeno di partenza, ma rimane comunque imprenscindibile la necessità di provare empiricamente, cioè sperimentando e lanciando il codice, quale sia effettivamente la migliore per il sistema e la GPU presa in considerazione.

Thread ID

Con le attuali nozioni, possiamo dunque realizzare il primissimo kernel:

__global__ void minimal( int* d_a)
{
    *d_a = 13;
}

Le sue uniche funzionalità sono però quelle di assegnare sempre alla stessa cella di memoria il valore 13, certamente uno spreco di potenza computazionale e di thread, considerando che questa operazione potrebbe facilmente venire richiamata da una configurazione di esecuzione con diverse centinaia di computazioni in parallelo.

Dunque, l'ultimo elemento base per poter realizzare un kernel effettivamente funzionante è capire come individuare l'identificatore (il numero) di un thread che sta operando su un determinato dato o insieme di dati. Per fare ciò, C for CUDA integra delle variabili sempre accessibili all'interno di funzioni di tipo __global__ e __device. Con sempre accessibili intendiamo il fatto che non c'è bisogno di alcuna dichiarazione per poter utilizzare queste variabili, che sono sempre attive e inizializzate ai valori corretti all'interno delle funzioni che stiamo scrivendo.

Esse sono:

• dim3 gridDim - ritorna le dimensioni della griglia in numero di blocchi. Poiché la grandezza della griglia è espressa in 2 dimensioni, il terzo valore è sempre 1.
• dim3 blockDim - ritorna le dimensioni del blocco in numero di thread.
• dim3 blockIdx - ritorna l'indice del blocco corrente.
• dim3 threadIdx - ritorna l'indice del thread corrente.

Attenzione alle maiuscole (trattandosi di C, è importante che la sintassi sia quella riportata). Si tratta in tutti i casi di variabili di tipo dim3, e dunque tridimensionali. Per accedere ai differenti campi, si aggiunge al nome della variabile l'identificatore del campo richiesto:

• blockDim.x: ritorna la dimensione sull'asse x del blocco (numero di thread guardando all'asse x).
• threadIdx.y: ritorna l'indice del thread corrente muovendosi sull'asse y.

Vogliamo annotare che con x si identifica l'asse orizzontale e con y l'asse verticale, esattamente come avviene sul piano cartesiano. Nel caso di spazi tridimensionali, ovviamente l'asse z è quello della profondità.

Per determinare l'indice di uno specifico thread all'interno di un kernel, la formuletta base, e che in generale viene sempre usata all'interno dei programmi CUDA è:

Tid = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.y * blockDim.x;

Per esempio, considerando i seguenti dati:

blockDim.x = 2;
blockDim.y = 4;
blockDim.z = 6;

threadIdx.x = 1;
threadIdx.y = 2;
threadIdx.z = 5;

L'indice di thread è:

Tid = 1 + 2 * 2 + 5 * 2 * 4 = 1 + 4 + 40 = 45

Finalmente, un kernel che esegue l'assegnazione di valori ad un vettore, tutte in parallelo tra loro, può essere scritto:

__global__ void assign( int* d_a, int value)
{
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    d_a[idx] = value;
}

Note sull'esecuzione dei kernel

Ci sono alcuni particolari da tenere in considerazione quando si stabilisce la configurazione di esecuzione di un kernel. Come detto, le griglie sono ordinate bidimensionalmente, mentre i blocchi possono essere tridimensionali. Per le griglie, trattandosi di blocchi, che non devono essere necessariamente attivi tutti allo stesso tempo, le dimensioni massime possono essere di 65536 elementi sia sull'asse X che sull'asse Y, con il massimo numero di blocchi che può raggiungere il prodotto tra i due. Per i blocchi, invece, si possono avere fino a 512 elementi sull'asse X, 512 sull'asse Y e 64 sull'asse Z, ma il prodotto delle combinazioni non può superare i 512 thread. Questo vale per tutte le versioni di hardware attualmente disponibili da parte di NVIDIA, ma non escludiamo che in futuro questi dati possano cambiare per nuovi prodotti. Diverso invece il discorso sui thread che possono essere attivi contemporaneamente su un singolo multiprocessore in un dato istante di tempo, ma ci inoltreremo in questo aspetto successivamente, parlando delle ottimizzazioni per massimizzare l'efficienza dei programmi CUDA.

Vogliamo anche sottolineare alcuni aspetti che normalmente non sono riportati nelle classiche guide di programmazione. In primo luogo, i kernel sono eseguiti sequenzialmente tra loro, ma una volta che la CPU lancia un kernel è libera di fare altre operazioni mentre esso viene eseguito dalla GPU. La chiamata è, cioè, asincrona, ed è dunque possibile far eseguire alla CPU altre parti di codice in parallelo alla GPU. Con le memcopy invece ciò ovviamente non accade poiché entrambi gli elementi di elaborazione sono coinvolti nell'operazione, che dunque è sincrona. Se si vuole mettere la CPU in attesa della terminazione di tutte le operazioni in esecuzione sul device, si può usare la primitiva:

cudaThreadSynchronize();

che attende la terminazione di tutte le chiamate CUDA precedenti.
L'asincronicità delle chiamate ai kernel rispetto alla CPU può portare anche a situazioni dove, se ci sono più programmi che stanno accedendo ad un singolo device, le invocazioni si accodano tra di loro. Visto che non c'è modo di sapere se il device sta attualmente eseguendo un kernel, questo potrebbe essere un piccolo problema nel caso si stia cercando di misurare i tempi di esecuzione. Motivo per il quale potrebbe essere consigliabile, quando si sta operando su un sistema di GPGPU computing condiviso per raccogliere dati prestazionali, usare adeguatamente le cudaThreadSynchronize(). C'è anche da tenere conto che il primo caricamento di un kernel comporta un certo overhead in termini di tempo che non è facilmente quantificabile e dipende molto anche dal sistema nel suo complesso. Consigliamo perciò di iterare numerose volte il kernel quando si vogliono registrare le prestazioni, per eliminare tali variabilità.  Anche l'occupazione del kernel in memoria non è facilmente intuibile, ma fortunatamente è di norma molto limitata visto che i kernel di buona qualità sono normalmente compatti. Infine, un ultima nota per quanto riguarda l'allocazione dei dati memoria. Le cudaMalloc() non dicono quanto spazio è rimasto disponibile, ma solo se l'allocazione è andata a buon fine o no. Dunque, se altre persone stanno usando lo stesso device (ad esempio, su un sistema di supercomputing), senza che ci sia noto, potremmo trovarci in situazioni nelle quali l'allocazione non va a buon fine pur avendo teoricamente a disposizione sulla scheda grafica ben più della memoria richiesta.

Il primo programma con un kernel funzionante

Mettiamo finalmente insieme tutti gli aspetti finora trattati, e scriviamo il nostro primo programma CUDA che compia realmente qualche operazione interessante.

// include di sistema
#include <stdio.h>
#include <assert.h>

// Semplice funzione per verificare gli errori in cuda
void checkCUDAError(const char *msg);

// implementazione del primo kernel
__global__ void primoKernel(int *d_a)
{
    // calcolo dell'indice di thread
    int idx = blockIdx.x*blockDim.x + threadIdx.x;
    
    // operazione: ad ogni elemento del vettore, assegna il valore dell'indice di blocco moltiplicato per mille più l'indice di thread
    d_a[idx] = 1000*blockIdx.x + threadIdx.x;  
}

// Dichiariamo il main
int main( int argc, char** argv) 
{
    // puntatore per la struttura dati sull'host
    int *h_a;

    // puntatore per la struttura dati sul device
    int *d_a;

    // definizione della grandezza della griglia e dei blocchi
    int numBlocks = 8;
    int numThreadsPerBlock = 8;

    // Allocazione della memoria sull'host
    size_t memSize = numBlocks * numThreadsPerBlock * sizeof(int);
    h_a = (int *) malloc(memSize);

    // e sul device
    cudaMalloc( (void **) &d_a, memSize );

    // Lancio del kernel
    dim3 dimGrid(numBlocks);
    dim3 dimBlock(numThreadsPerBlock);
    primoKernel<<< dimGrid, dimBlock >>>( d_a );

    // blocca la CPU fino al completamento del kernel sul device
    cudaThreadSynchronize();

    // controlla se l'esecuzione del kernel ha generato qualche errore
    checkCUDAError("kernel execution");

    // Esegue la copia dei risultati dalla memoria del device a quella dell'host
    cudaMemcpy( h_a, d_a, memSize, cudaMemcpyDeviceToHost );

    // controlla se l'esecuzione della memcopy ha generato qualche errore
    checkCUDAError("cudaMemcpy");

    // Verifica che i dati ritornati all'host siano corretti
    for (int i = 0; i < numBlocks; i++)
    {
        for (int j = 0; j < numThreadsPerBlock; j++)
        {
            assert(h_a[i * numThreadsPerBlock + j] == 1000 * i + j);
        }
    }

    // libera la memoria sul device
    cudaFree(d_a);

    // libera la memoria sull'host
    free(h_a);

   // Se il programma arriva fin qua, tutto corretto!
    printf("Corretto!\n");

    return 0;
}

void checkCUDAError(const char *msg)
{
    cudaError_t err = cudaGetLastError();
    if( cudaSuccess != err) 
    {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString( err) );
        exit(-1);
    }                         
}