Nella precedente parte di questo tutorial, abbiamo scritto il nostro primo kernel in C per CUDA. In questa terza puntata, completeremo le ultime basi di programmazione, che ci permetteranno di essere operativi in maniera abbastanza efficiente con l'ambiente NVIDIA. In particolare, vedremo come gestire l'allocazione di memoria in 2D e come poter sfruttare la shared memory. Faremo anche un brevissimo excursus sui tipi vettoriali integrati nell'ambiente, le istruzioni di sincronizzazione atomiche supportate, la gestione degli errori e il funzionamento di NVCC. Conclusa questa parte, sarà possibile poi passare alla parte di ottimizzazione, che permetterà effettivamente di ottenere le prestazioni dal proprio codice CUDA.
Gestione della memoria 2D
È piuttosto evidente che il modello di programmazione di CUDA si adatti molto bene ad algoritmi che operino su matrici. L'organizzazione dei thread in blocchi (tridimensionali) e dei blocchi in griglie (bidimensionali) ne è un chiara dimostrazione. In questa situazione, è sicuramente molto gradito avere funzioni apposite per gestire adeguatamente la memoria come uno spazio bidimensionale. Del resto, gestire una matrice come un vettore, con la conseguente necessità di linearizzarla, spesso non è neanche un'operazione così naturale e facile.
NVIDIA mette a disposizione in CUDA alcune primitive per allocare la memoria in 2D, che sono anche strettamente legate alla modalità operativa dell'hardware, e ne favoriscono in un certo modo le prestazioni. In particolare, in questa parte del nostro tutorial vogliamo focalizzare l'attenzione sulla funzione:
cudaMallocPitch(void** devPtr, size_t* pitch, size_t widthInBytes, size_t height);
Tramite questa funzione, possiamo allocare in memoria una struttura, puntata da devPtr, che ha righe della dimensione widthInBytes (come il nome fa intuire, la grandezza va specificata in byte) ed è alta un numero di linee equivalente a height. Oltre a ritornare l'indirizzo di partenza di della struttura in in devPtr, questa primitiva torna anche un valore, in bytes, di pitch. Il pitch indica la grandezza che dovrebbe avere una linea per ridurre, il più possibile, il numero di transazioni di memoria eseguite dall'hardware per accedere ai dati stessi.
Esiste anche una apposita funzione di memcopy:
cudaMemcpy2D(void *dst, size_t dpitch, const void *src, size_t spitch, size_t width, size_t height, enum cudaMemcpyKind direction);
Che prende copia, sulla destinazione *dst, una struttura *src, rispettivamente con pitch dpitch e spitch, con linee di lunghezza width su height linee.
Se ancora non vi è chiaro il significato del pitch, la seguente figura dovrebbe aiutare a chiarirsi le idee:
Supponiamo (come effettivamente è, vedremo meglio nella parte di ottimizzazioni quando parleremo molto più dettagliatamente dell'architettura delle attuali GeForce) che in una singola transazione di memoria sia possibile caricare 16 elementi da 32 bit (64 byte). Se la nostra matrice è formata da linee di solo da 15 elementi, con una normale allocazione esse saranno poste tutte una di seguito all'altra. Prendendo il caso, più semplice, delle schede GeForce 8 e 9 con compute capability 1.0 e 1.1, questo significa che il controller di memoria eseguirà almeno:
- 1 transazione da 64 byte per recuperare la prima linea
- 1 transazione da 64 byte per recuperare il primo valore della seconda linea
- 1 transazione da 64 byte per recuperare i restanti 14 valori della terza linea
Il disallineamento poi continuerà a protrarsi per tutte le successive linee della matrice.
Se, invece, si fa in modo di aggiungere un elemento “vuoto” a ciascuna linea della matrice, il controller eseguirà:
- 1 transazione da 64 byte per recuperare la prima linea
- 1 transazione da 64 byte per recuperare la seconda linea
E via dicendo per tutto il resto della matrice.
Il pitch esprime giustappunto in byte la lunghezza che dovrebbe avere la linea della matrice per raggiungere l'allineamento di memoria che riduca il più possibile il numero di transazioni per il caricamento dei dati. In questo esempio, abbiamo linee da 60 byte (15 elementi per 32 bit), dunque il pitch sarà 64 byte. Nel caso la nostra linea sia, per esempio, 116 byte (29 elementi da 32 bit), allora il pitch sarà di 128 byte, e così via.
Con la cudaMallocPitch la memoria viene proprio allocata con le righe della lunghezza “giusta” per ridurre il numero di transazioni. Questo significa, che nella matrice allocata nella memoria del device, la nuova linea inizierà dopo “pitch” byte. Ciò influisce direttamente nel modo in cui andremo ad accedere alla matrice nel nostro kernel. Infatti, un kernel minimale 2D che usa matrice:
__global__ void assign2D(int* d_a, int row, int value)
{
int iy= blockDim.y * blockIdx.y + threadIdx.y;
int ix= blockDim.x * blockIdx.x + threadIdx.x;
int idx = iy * row + ix;
d_a[idx] = value;
}
…
assign2D<<<dim3(64,64), dim3(16,16)>>>(d_a, row, value);
Dove row rappresenta il numero di elementi su una singola riga, diventa:
__global__ void assign2D(int* d_a, int mat_pitch, int value)
{
int iy= blockDim.y * blockIdx.y + threadIdx.y;
int ix= blockDim.x * blockIdx.x + threadIdx.x;
int idx = iy * mat_pitch + ix;
d_a[idx] = value;
}
…
int mat_pitch=pitch/sizeof(int);
assign2D<<<dim3(64,64), dim3(16,16)>>>(d_a, mat_pitch,value);
Con il numero di elementi sulla singola riga rappresentata dal pitch ottenuto dalla primitiva cudaMallocPitch, diviso per la grandezza in byte del singolo elemento della struttura dati, in questo caso composta da interi. L'operazione di divisione non va dimenticata, visto che il valore di pitch ritornato è in byte e il calcolo dei puntatori all'interno del kernel non funzionerebbe più correttamente senza.
Qualificatori di variabili
Nella scorsa puntata, abbiamo visto i qualificatori per le funzioni. In CUDA esistono anche dei qualificatori di variabile, che permettono sostanzialmente di specificare in quale tipo di memoria devono essere allocate le variabili dichiarate in un kernel.
In particolare, possiamo distinguere i qualificatori:
- __device__
Salva nella device memory (lenta, senza cache)
Allocata con cudaMalloc
Accessibile da tutti i thread
ha come tempo di vita quello dell'applicazione
- __shared__
salvato nella shared memory on chip (veloce)
Allocata durante la configurazione di esecuzione o durante la compilazione
Accessibile da tutti i thread nello stesso blocco
ha come tempo di vita quello dell'esecuzione del kernel
E' possibile anche non qualificare le variabili. Di base, se ad una dichirazione di variabile non si aggiunge nulla, il compilatore agisce:
- salvando gli scalari e i le strutture vettoriali predefinite nei registri
- salvando array di più di quattro elementi nella memoria del device
Tipi vettore predefiniti
I tipi di vettore predefiniti in CUDA possono essere usati sia nel codice per la GPU che nel codice per la CPU. Uno di questi tipi è naturalmente il dim3, che ha tre campi e permette di specificare le dimensioni ed ha un valore di default di (1,1,1).
Più in generale, abbiamo tipi vettoriali predefiniti per char, short, integer, long e float. dim3 è basato su uint3. I primi tre hanno le varietà senza segno:
[u]char[1...4], [u]short[1...4], [u]int[1...4], [u]long[1...4], float[1...4];
Si accede alle strutture con i campi x, y, z e w:
uint4 param;
int y = param.y;
Usare la shared memory
Per usare la shared memory si deve ovviamente qualificare la variabile che si vuole allocare in tale memoria con la particella __shared__. Poiché questa memoria è particolarmente veloce, il suo uso è consigliato per allocare dei vettori nei quali copiare quei dati che i thread di un blocco:
- accedono con grande frequenza
- usano per operazioni complesse
Un array in shared memory è visibile da tutti (e soli) i thread all'interno di un singolo blocco. La shared memory è limitata a 16 KB per ogni multiprocessore. Questo significa che in un dato istante, su un multiprocessore dell'architettura potranno essere allocati al massimo tanti blocchi quanti permessi dalla shared memory. Come per il discorso sul numero di thread attivi, ciò non significa che non è possibile avere nel kernel più blocchi di questa cifra. Semplicemente, in un dato momento, i thread che il multiprocessore potrà eseguire saranno scelti tra quelli attivi. Gli altri verranno assegnati solo in momenti successivi, quando quelli precedenti saranno terminati. La shared memory non è ovviamente l'unico vincolo per determinare il massimo numero di blocchi (e thread) attivi su un multiprocessore in un certo momento: gli altri due fattori sono il massimo numero di thread gestibili dallo scheduler e il numero di registri (per thread) che possono essere usati. Vedremo meglio nella parte di ottimizzazione come agire per cercare di identificare le migliori configurazioni per massimizzare l'efficienza.
Utilizzare la shared memory significa poter accedere ai dati salvati in essa con prestazioni e latenze praticamente equivalenti a quelle dei registri, dunque elevatissime. Nelle parte di ottimizzazioni vedremo però che l'organizzazione a banchi di tale memoria consegue nella necessità di accedere ad essa secondo schemi ben precisi per poter ottenere le migliori performance, eliminando i conflitti tra i banchi. Il vantaggio che la shared memory è in grado di dare può, comunque, in molti casi essere quantificabile con un intero ordine di grandezza, in quanto la sua funzione è quella, principalmente, di aumentare considerevolmente la banda a disposizione per i thread.
Esistono due metodi per allocare shared memory:
se la grandezza delle strutture dati da allocare è conosciuta al momento della compilazione:
__global__ void kernel(...)
{
…
__shared__ float sData[256];
…
}
int main(void)
{
…
kernel<<<nBlocks, blockSize>>>(...);
…
}
se invece è conosciuta al momento del lancio del kernel:
__global__ void kernel(...)
{
…
extern __shared__ float sData[];
…
}
int main(void)
{
…
int smBytes = blockSize*sizeof(float);
kernel<<<nBlocks, blockSize, smBytes>>>(...);
…
}
Viene cioè aggiunto un terzo parametro alla configurazione di esecuzione. La modalità di allocazione con grandezza impostata al lancio del kernel ha però delle limitazioni. È infatti possibile dichiarare una sola struttura dati, la cui grandezza è definita nell'unico parametro di configurazione a disposizione, e dunque diventa necessario gestirsi manualmente la sua organizzazione se in essa si vogliono depositare dati di tipologie diverse. Nel caso della grandezza conosciuta a compile time, invece, è possibile dichiarare tutte le strutture volute.
Sincronizzazione e operazioni atomiche
Ovviamente, scrivere in una memoria condivisa può portare a condizioni di corse critiche: nulla vieta di leggere con dei thread dei dati dalla memoria globale, scriverli nella shared memory, e poi riutilizzare gli stessi thread, magari con uno schema d'accesso diverso dal precedente, per operare su quei dati che si sono copiati nella shared memory. Anzi, questa sarà, con la massima probabilità, la cosa comunemente fatta utilizzando la shared memory. Si rivela dunque necessaria una primitiva che sia in grado di bloccare i thread finchè tutti non la raggiungono, in modo da assicurarsi che tutti abbiano, per esempio, terminato le operazioni di scrittura. Tale funzione, che nell'ambito della programmazione multithreaded è normalmente conosciuta come barriera, in CUDA è realizzata tramite la primitiva:
__syncthreads();
Questo costrutto sincronizza tutti i thread all'interno di un blocco, cioè tutti quei thread che, appunto, possono accedere alle stesse locazioni in shared memory. D'altra parte, thread di blocchi differenti non vengono sincronizzati tra loro, così come non possono collaborare tramite shared memory. Quando viene eseguito tale costrutto, nessun thread del blocco può superarlo finchè anche tutti gli altri non lo hanno raggiunto. Il suo utilizzo è principalmente quello di impedire conflitti (lettura dopo scrittura - RAW, scrittura dopo lettura - WAR, scrittura dopo scrittura – WAW) quando si accede alla shared memory. Da notare che il suo utilizzo in codice condizionale (for, if) è permesso solo se la condizione è uniforme per tutto il blocco.
CUDA permette anche l'utilizzo di alcune operazioni atomiche, operazioni, cioè, che vengono eseguite da un thread alla volta perchè se eseguite in parallelo comporterebbero dei problemi di sincronizzazione. L'esempio classico è quello dall'addizione atomica, dove si vuole che ciascun thread aggiunga un valore ad un valore precedente, e ovviamente l'operazione deve essere fatta sequenzialmente per permettere a un thread di leggere l'incremento fatto precedemente da un altro thread. Bisogna fare attenzione al fatto che il supporto alle istruzioni atomiche è stato introdotto solo sulle schede con compute capability 1.1 (sostanzialmente, i chip usciti dopo le prime GeForce 8800 GTX e GT, che erano basate su G80: G84, G92, ecc).
In dettaglio, è possibile eseguire operazioni atomiche di tipo:
- associativo, su interi con segno e senza
- addizione, sottrazione, minimo, massimo
- and, or, xor
- incremento, decremento
- Sosistituzione (exchange), comparazione, e scambio (swap)
Bisogna anche sottolineare che con i chip compute capability 1.1 queste funzioni sono supportate solo in global memory (quindi piuttosto lente) e su interi a 32 bit, mentre con la compute capability 1.2 sono supportate anche in shared memory, quindi a livello di blocco, e su interi a 64 bit. Bisogna precisare che non esistono attualmente chip con compute capability 1.2, ma le schede basate su architettura GT200, cioè tutte le GeForce della serie GTX 2xx, sono compute capability 1.3, che supporta ovviamente tutte le funzionalità 1.2 e aggiunge la doppia precisione in virgola mobile.
Riportare gli errori alla CPU
Tutte le chiamate CUDA, fatta eccezione per i lanci di kernel, ritornano un codice di errore del tipo:
cudaError_t type;
Esiste una funzione che può essere sfruttate nel codice per riportare l'errore all'utente:
cudaError_t cudaGetLastError(void);
ritorna il codice dell'ultimo errore (anche il nessun errore ha un codice), e può essere usato per ottenere gli errori nell'esecuzione di un kernel. Combinandola con:
char* cudaGetErrorString(cudaError_t code);
che ritorna una stringa di caratteri che descrive l'errore, si può fornire all'utente una descrizione più significativa del problema:
printf(“%s\n”, cudaGetErrorString(cudaGetLastError()));
Un kernel con shared memory
Riportiamo qua l'esempio di un kernel che usa blocchi e shared memory. Lo scopo del kernel è quello di leggere un array e scriverlo invertito.
// include di sistema
#include <stdio.h>
#include <assert.h>
// Semplice funzione per verificare gli errori in cuda
void checkCUDAError(const char* msg);
// implementazione del kernel
__global__ void reverseArrayBlock(int *d_out, int *d_in)
{
// dichiarazione del vettore da allocare in shared memory
// allocato con determinazione delle dimensioni al lancio del kernel
extern __shared__ int s_data[];
// calcolo dell'offset di blocco di input
int inOffset = blockDim.x * blockIdx.x;
// calcolo dell'indice dell'elemento del vettore letto da un singolo thread
int in = inOffset + threadIdx.x;
// Legge un elemento per ciascun thread dalla memoria del device
// lo salve in ordine inverso nella shared memory
s_data[blockDim.x - 1 - threadIdx.x] = d_in[in];
// Blocca i thread del blocco
// finchè non tutti hanno terminato l'operazione di scrittura
__syncthreads();
// scrive i dati dalla memoria shared in ordine diretto,
// ma usando l'ordine inverso per i blocchi:
// il primo blocco diventa l'ultimo e così via
int outOffset = blockDim.x * (gridDim.x - 1 - blockIdx.x);
int out = outOffset + threadIdx.x;
d_out[out] = s_data[threadIdx.x];
}
// dichiarazione del main
int main( int argc, char** argv)
{
// puntatore alla memoria host e grandezza
int *h_a;
int dimA = 256 * 1024; // 256K elementi (1MB)
// puntatore alla memoria del device
int *d_b, *d_a;
// grandezza del blocco in thread
int numThreadsPerBlock = 256;
// Calcolo del numero di blocchi in base al numero di thread per blocco
int numBlocks = dimA / numThreadsPerBlock;
// calcolo della quantità, in byte, di shared memory richiesta per blocco
int sharedMemSize = numThreadsPerBlock * sizeof(int);
// allocazione della memoria sull'host e sul device
size_t memSize = numBlocks * numThreadsPerBlock * sizeof(int);
h_a = (int *) malloc(memSize);
cudaMalloc( (void **) &d_a, memSize );
cudaMalloc( (void **) &d_b, memSize );
// Inizializza l'array iniziale sull'host
for (int i = 0; i < dimA; ++i)
{
h_a[i] = i;
}
// Copia l'array dall'host al device
cudaMemcpy( d_a, h_a, memSize, cudaMemcpyHostToDevice );
// lancio del kernel
dim3 dimGrid(numBlocks);
dim3 dimBlock(numThreadsPerBlock);
reverseArrayBlock<<< dimGrid, dimBlock, sharedMemSize >>>( d_b, d_a );
// blocca fino al completamento del device
cudaThreadSynchronize();
// Controlla se l'esecuzione del kernel ha generato un errore
checkCUDAError("kernel invocation");
// copia dei risultati dal device all'host
cudaMemcpy( h_a, d_b, memSize, cudaMemcpyDeviceToHost );
// Controllo errori CUDA
checkCUDAError("memcpy");
// verifica che i dati ritornati all'host siano corretti
for (int i = 0; i < dimA; i++)
{
assert(h_a[i] == dimA - 1 - i );
}
// libera la memoria sul device
cudaFree(d_a);
cudaFree(d_b);
// 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(EXIT_FAILURE);
}
}
Compilare le applicazioni C per CUDA
Le applicazioni C per CUDA possono essere composte da file con codice C/C++ standard e file con codice che sfrutta le estensioni CUDA (tutti i file .cu che abbiamo usato negli esempi). I file con solo codice C/C++ possono essere dati in pasto ad un compilatore standard, i file invece che contengono anche estensioni CUDA devono passare sotto NVCC. In entrambi i casi, si generano dei file oggetto che poi il linker integra, permettendo di ottenere un singolo eseguibile con i binari per CPU e GPU.
Più in dettaglio, anche nei file .cu è contenuto codice C per la CPU e codice C per CUDA per la GPU. NVCC, che è un compiler driver, separa i due codici, lanciando il compilatore di sistema (Visual C per Windows, o GCC per Linux/Unix) per la parte per la CPU e invocando il CUDA compiler per la parte destinata alla GPU. Il CUDA compiler genera un binario di tipo PTX (Parallel Thread eXecution), che rappresenta l'Instruction Set Virtuale per i chip grafici NVIDIA, definendone anche il modello di programmazione, le risorse di esecuzione e lo stato. Un ulteriore traslatore (o il compilatore stesso, o l'interprete a runtime di CUDA) trasforma il codice PTX nel codice binario dell'architettura fisica target.
Tutti gli eseguibili con codice CUDA richiedono due librerie dinamiche:
- La CUDA core library (cuda)
- La CUDA runtime library (cudart)
Quest'ultima, nel caso sia usata l'API a runtime (che permette di caricare codice PTX che viene compilato a runtime), carica autonomamente la core library.
Debugging con la modalità di emulazione
Come abbiamo precedentemente visto, un eseguibile compilato in modalità di emulazione (nvcc -deviceemu), funziona sull'host usando il runtime di CUDA. In tale situazione, un singolo thread è emulato da un thread dell'host, dunque dal punto di vista delle prestazioni è impossibile avere indicazioni significative. D'altra parte, ciò significa anche che non c'è bisogno di avere alcun driver CUDA e alcun chip che supporti CUDA per poter scrivere e provare il codice. E si possono utilizzare modalità molto classiche per la ricerca degli errori. Innanzitutto, si possono usare i programmi standard di debugging sull'host (ad esempio gdb, che permette di impostare breakpoint). È poi possibile ispezionare i dati diretti al device direttamente dall'host, e usare funzioni disponibili sull'host dal codice del device. Pensiamo, ad esempio, alle funzioni di input/output come le printf. La sequenzializzazione può anche permettere di rilevare delle corse critiche causate da un uso sbagliato delle _syncthreads(). D'altra parte, la stessa sequenzializzazione può portare a risultati differenti nel caso ci siano accessi simultanei di più thread ad una stessa locazione di memoria. Inoltre, è si possibile accedere a puntatori sul device dall'host, o viceversa, in modalità di emulazione, ma tali operazioni non possono più essere attuate quando il codice per la GPU viene eseguito sul device.
