In occasione del Pi day (il 14 marzo), il mio precedente post sul calcolo del pi greco col metodo Monte Carlo ha ricevuto nuova attenzione. In particolare un utente incontrato su un gruppo facebook molto popolare su Arduino, @carmelopellegrino, ha voluto dare un importante contributo al codice multi thread, che nel frattempo avevo messo su github. Nella seconda parte invece sono presentati i risultati di una sperimentazione con CUDA, avvenuta grazie al prezioso contributo di @akiross.
La bellezza della condivisione, no?
Ma partiamo dall’inizio. Per prima cosa mi sono spostato su linux, per quanto windows sia bello e simpatico nel far sentire un utente a casa propria, linux è molto più adatto alla programmazione, c’è poco da fare. è un sistema operativo creato dai programmatori per i programmatori (il che porta ad un altro grosso difetto, ma magari ne parleremo prossimamente). Detto questo, il passaggio si è portato dietro molti più problemi di quanti pensassi.
Prima di tutto, per la versione multithread avevo scelto di usare il sistema dei pthread, cioè i thread di POSIX. Non voglio scendere nel dettaglio di cosa sia POSIX e cosa siano i pthread perché non sono la persona adatta a farlo. In due parole: POSIX è uno standard che principalmente definisce diversi aspetti dei sistemi operativi Unix-like. Tra questi c’è anche la gestione dei thread, e i thread che seguono le regole di POSIX prendono il nome di pthread. Tra i sistemi operativi Unix-like c’è anche Linux, che è il motivo per cui stiamo affrontando questo discorso. Possiamo capire quindi che i pthread funzionano molto più nativamente su Linux rispetto a Windows, ragione per cui ero portato a credere che avrei guadagnato in performance. In parte era vero, ma solo dopo alcune importanti modifiche al codice.
La prima modifica, che non influenza le prestazioni, riguarda essenzialmente la misura del tempo di esecuzione. Il metodo di misura che usavo su windows era:
clock_t time = clock();
// codice di cui misurare il tempo...
//...
time = clock() - time;
printf("Total CPU time: %lf seconds\n", (double)time/(CLOCKS_PER_SEC));
Essenzialmente si salva il valore di clock() prima dell’esecuzione della parte di cui misurare il tempo, finita l’esecuzione si ricava la differenza di clock() tra inizio e fine, e si divide per CLOCKS_PER_SEC, ovvero una costante che contiene il numero di clock che vengono eseguiti per secondo nel nostro sistema. Su windows questo codice funziona senza problemi, ma linux nasconde alcune insidie su quest’aspetto.
Nonostante clock() sia una funzione standard del C, sui due sistemi operativi misura due cose diverse: su windows misura il clock time (ovvero il tempo d’orologio), su linux invece misura il cpu time (ovvero il tempo di esecuzione sulla cpu). Questo normalmente non comporta particolari differenze, se non quando si esegue un programma su più thread e si vuole ottenere il tempo dell’orologio trascorso durante l’esecuzione. Su linux infatti ogni thread può occupare il 100% della cpu, portandoci ad un uso totale dell’800% nel caso di 8 thread che occupano ciascuno il 100% (sembra un nonsense, come posso usare più cpu di quella che ho?). Comunque si tratta di un problema risolvibile con estrema facilità: basta dividere il tempo totale per il numero di thread, e otterremo il nostro tempo d’orologio. Nonostante la facilità con cui può essere risolto, ha richiesto qualche tempo per capire cosa ci fosse dietro alle misurazioni di tempo sbagliate, per cui ho pensato che valesse la pena di sottolineare quest’aspetto.
L’altro piccolo problema riguardava la funzione Beep(), che generava il tono a computazione terminata. Questa funzione è contenuta nella libreria windows.h, che ovviamente non è disponibile su linux. Per ovviare a problemi di portabilità, ho deciso di utilizzare il carattere speciale \a, che prende il nome di “bell character”. Questo carattere, quando stampato con una printf() produce un suono che attira l’attenzione dell’utente. Il nome peraltro è curioso: è chiamato così perché nelle macchine per corrispondenza telematica (le telescriventi) azionava un vero e proprio campanello elettromeccanico.
Purtroppo questo non era tutto, c’erano dei bug che minavano pesantemente la velocità di esecuzione, che Carmelo ha scovato e mi ha aiutato a correggere.
Il principale problema si insidiava nel fatto che la funzione rand() non è thread-safe. Il che significa, essenzialmente, che durante la sua esecuzione esegue delle operazioni su una porzione di memoria condivisa. Il problema è che a questa memoria condivisa può accedere solamente un thread per volta, per non causare collisioni o errori. Questo comportava un “accodamento” dei thread per accedere alla memoria, il che significava che molti thread per molto tempo rimanevano bloccati in attesa di accedere alla memoria. Se non gestita adeguatamente dal sistema operativo poi, questa situazione può probabilmente causare problemi ancora peggiori.
La funzione rand() inoltre non è rientrante, il che significa che il sistema non può interromperla e poi riprendere da dove era rimasto (appunto perché nel frattempo un altro thread potrebbe accedere e modificare quella regione di memoria).
Per rendere questo problema più chiaro, facciamo un esempio: supponiamo che la funzione rand() possa avere quattro stati interni: A, B, C e D. Supponiamo che si “ricordi” in che stato è, registrandolo in una variabile contenuta nella porzione di memoria condivisa. Supponiamo inoltre che proceda ordinatamente attraverso questi quattro stati: A -> B -> C -> D. Quando il thread1 chiama questa funzione, lei si dispone dapprima nello stato A, attraversato lo stato A passa allo stato B, e registra l’attuale stato in memoria. Supponiamo ora che thread1 venga interrotto da un interrupt di sistema (ad esempio generato dal fatto che l’utente muove il mouse ed è richiesta una certa computazione). Non avendo thread liberi a disposizione, il sistema, interrompe l’esecuzione di rand() di thread1. A questo punto, mentre thread1 esegue la richiesta del sistema, subentra la chiamata a rand() di thread2, che precedentemente era in attesa di accedere alla memoria condivisa. Questa entra nello stato A, e viene scritto nella porzione di memoria condivisa. A questo punto, sia che thread1 riprenda immediatamente l’esecuzione, sia che thread2 arrivi al termine, quando thread1 riprenderà controllo della memoria, il suo stato e possibilmente altre variabili avranno valori inaspettati, questo può generare problemi.
Questo esempio potrebbe essere il caso della rand(), oppure no. Sta di fatto che, come è stato ampiamente menzionato su internet da persone ben più esperte di me ([1], [2]), la funzione rand() soffre di questi problemi, pertanto non è corretto usarla in più thread simultaneamente.
Quindi cosa usare al suo posto? Carmelo ha suggerito di usare la funzione rand_r() che è appositamente studiata per essere rientrante e thread-safe.
Questo ha portato a una nuova versione del codice:
typedef struct _thread_data_t {
int tid; // thread id given by the order of creation
int belongs; // stores here the number of points that belong the circle
unsigned int reentrantSeed; //seed for the reentrant random function rand_r()
} thread_data_t;
void *thr_func(void *arg) {
thread_data_t *data = (thread_data_t *)arg;
// making sure that every thread gets a different set of random numbers
data->reentrantSeed=time(NULL)+data->tid;
printf("Thread %d is running\n", data->tid); // just to say hi
static const unsigned long int randmax_big=(((RAND_MAX+1.0)*(RAND_MAX+1.0))-1.0);
double x,y,x1,x2,y1,y2;
unsigned int executed=0; // stores the number of executed points
data->belongs=0; //initalizes to 0 the counter of matched points
// go through each point
for(executed;executed<THREAD_POINTS;executed++){
x1=rand_r(&data->reentrantSeed);
x2=rand_r(&data->reentrantSeed);
y1=rand_r(&data->reentrantSeed);
y2=rand_r(&data->reentrantSeed);
//the enhanced random formula
x=x1*RAND_MAX_P1+x2;
x=x/randmax_big;
y=y1*RAND_MAX_P1+y2;
y=y/randmax_big;
//chech if it belongs to the circle or not
if(x*x+y*y<1){
data->belongs++;
}
}
// say goodbye. the datas are stored in the thr_data array
pthread_exit(NULL);
}
A differenza della rand(), la rand_r() richiede che come argomento gli sia passato il seed del generatore. Il seed è quindi conservato nella struct thread_data_t. Un array di queste struct (una per ogni thread) è creato dal main per assegnarvi i thread ID, e successivamente passata per indirizzo al thread. Questa implementazione è risultata comunque problematica. Per comprendere esattamente il perché, dobbiamo scavare ancora un po’ più a fondo nel funzionamento dei processori.
I più attenti avranno notato, magari scorrendo il sito internet del rivenditore di fiducia, che quando ci si addentra nelle specifiche dei processori, molto spesso vengono menzionate delle memorie cache, a volte accompagnate da diciture come L1, L2 o L3, espresse in megabyte. Per capire i successivi problemi, è necessario avere una chiara idea di cosa sono e come funzionano le memorie cache dei processori.
Partiamo dal che cosa è una memoria cache. Come avrete ormai capito non sono un maestro delle definizioni, me la cavo meglio a spiegare come funzionano le cose, ma provo a darvi comunque una definizione. Una cache è una piccola memoria, straordinariamente veloce (molto più della ram) che viene usata per conservare dati spesso richiesti dal processore. Per velocizzarne la lettura vengono conservati in queste memorie dati e istruzioni del programma. Ma non tutte le cache sono uguali.
Nei processori attuali troviamo ben 3 cache: L1, L2 e L3. La L che precede il numero sta per level, e il numero indica il vero e proprio livello della cache. Essenzialmente i livelli funzionano così: più il livello è alto (livello 1) e più la cache è veloce e piccola. Inoltre la cache L3 è condivisa tra i vari core, mentre degli altri due livelli ogni core ha la sua copia.
Il compito di amministrare la cache comunque, decidere quali dati ci vanno registrati e quali possono essere eliminati, è un compito estremamente delicato. È inoltre necessario sapere che nella cache non viene copiata solo la variabile richiesta, ma un intero blocco di memoria, che prende il nome di linea. Questo blocco di memoria contiene, oltre alla variabile richiesta, le variabili contigue. Spesso infatti è ragionevole ritenere che, se il processore necessita di accedere a una variabile, presto necessiterà di accedere anche a quelle vicine. Inoltre i dati modificati da un core vengono conservati per un certo periodo di tempo nella cache per consentire di fare successive modifiche, prima di copiarli di nuovo nella memoria RAM. Ciascun core può tenere traccia delle modifiche di un segmento di ram tramite un flag che viene chiamato dirty bit. Una linea non aggiornata nella ram viene infatti detta linea sporca. Il fatto che una linea sia sporca e che gli altri core possiedano una copia della stessa linea, implica la necessità di aggiornarla nella ram, da cui gli altri core avranno la possibilità di aggiornare la loro copia. Questo naturalmente implica un enorme dispendio di risorse computazionali e di tempo.
Perché abbiamo parlato di questo? Perché quasi per definizione, gli array occupano blocchi contigui di memoria. E come abbiamo detto, le struct che conservano i dati sono riunite in un array creato nel main. Ma se continuiamo a leggere e ad aggiornare i dati della struct, anche chi ha salvato dei dati contigui avrà necessità di aggiornarli in seguito alle modifiche di un solo thread, portando a un enorme spreco di risorse computazionali.
Un’ultima miglioria rispetto al codice originale, che sembra irrilevante ma è in grado di dare un contributo sensibile alle performance è questa: per ottenere una maggior quantità di numeri random, è stata usata questa formula: x = (x1 * RAND_MAX_P1 + x2) / randmax_big dove x1 e x2 sono due numeri casuali “piccoli”, generati da rand_r(), RAND_MAX_P1 è il massimo valore che possono assumere x1 e x2 e randmax_big è il massimo valore che può assumere il risultato della forumla. Dividiamo tutto per randmax_big per ottenere valori compresi tra zero e uno, da usare per le prove. Ebbene, sostituendo randmax_big con il suo inverso (calcolato una volta sola) e moltiplicando il risultato per l’inverso invece che dividere, otteniamo un miglioramento delle prestazioni. Matematicamente il risultato è lo stesso, ma a livello di computazione la divisione è più costosa.
Tutte queste migliorie hanno consentito un ottimo incremento delle performance dell’algoritmo, ma volevo spingermi più in là.
Già quando avevo ottenuto i primi risultati, prima della pubblicazione del primo articolo, con qualche compagno di corso si era parlato del fatto che sarebbe stato un ottimo codice da eseguire su una GPU. Perché su una GPU? Le CPU sono in grado di eseguire calcoli a velocità impressionanti, ma in modalità sequenziale. Già da qualche anno quasi tutte le CPU dispongono di più core, che consentono di eseguire simultaneamente processi diversi su un’unica CPU. Il compito delle GPU invece è radicalmente diverso. Nonostante entrambe le unità di elaborazione eseguano calcoli, nelle GPU bisogna eseguire ripetitivamente le stesse operazioni, relativamente semplici, per miliardi di volte, cioè per ciascun pixel. Questo ha portato le GPU ad evolversi nel campo della computazione parallela. Una GPU moderna infatti è in grado di eseguire decine di migliaia di thread contemporaneamente, a una condizione: che ogni thread esegua le stesse operazioni degli altri.
Dopo aver consultato un amico esperto e aver ricevuto qualche dritta da lui, mi sono convinto che la computazione su GPU avrebbe potuto fare al caso mio. In particolare il framework sviluppato da NVIDIA per le sue GPU, CUDA, mi è parso una buona soluzione. Poi in realtà tra l’università ed altri progetti in cui sono stato impegnato, non ho mai avuto occasione di mettermi a studiare veramente CUDA e la computazione su GPU. Fino a quando, qualche settimana fa, in una chat che frequento ho conosciuto Alessandro (@akiross). Lui CUDA lo usa per lavoro. Ho pensato di parlargli un po’ del mio problema (uno dei – cit.), per vedere se avrebbe saputo darmi qualche consiglio. Evidentemente il mio passatempo l’ha in qualche modo divertito, perché praticamente mi ha detto (non sono testuali parole): “sai che c’è? Dammi mezz’oretta che scrivo la mia implementazione, poi ne discutiamo“. Effettivamente qualche decina di minuti e una teiera bruciata dopo, se n’è uscito con questo gist. (il gist è uno strumento di github che permette di condividere bozze di codice non abbastanza significative da giustificare un repository tutto loro). Riporto il codice qui per comodità:
// Approximation of Pi using a simple, and not optimized, CUDA program
// Copyleft Alessandro Re
//
// GCC 6.x not supported by CUDA 8, I used compat version
//
// nvcc -std=c++11 -ccbin=gcc5 pigreco.cu -c
// g++5 pigreco.o -lcudart -L/usr/local/cuda/lib64 -o pigreco
//
// This code is basically equivalent to the following Python code:
//
// def pigreco(NUM):
// from random import random as rand
// def sqrad():
// x, y = rand(), rand()
// return x*x + y*y
// return 4 * sum(1 - int(test()) for _ in range(NUM)) / NUM
//
// Python version takes, on this machine, 3.5 seconds to compute 10M tests
// CUDA version takes, on this machine, 1.6 seconds to compute 20.48G tests
//
#include <iostream>
#include <limits>
#include <cuda.h>
#include <curand_kernel.h>
using std::cout;
using std::endl;
typedef unsigned long long Count;
typedef std::numeric_limits<double> DblLim;
const Count WARP_SIZE = 32; // Warp size
const Count NBLOCKS = 640; // Number of total cuda cores on my GPU
const Count ITERATIONS = 1000000; // Number of points to generate (each thread)
// This kernel is
__global__ void picount(Count *totals) {
// Define some shared memory: all threads in this block
__shared__ Count counter[WARP_SIZE];
// Unique ID of the thread
int tid = threadIdx.x + blockIdx.x * blockDim.x;
// Initialize RNG
curandState_t rng;
curand_init(clock64(), tid, 0, &rng);
// Initialize the counter
counter[threadIdx.x] = 0;
// Computation loop
for (int i = 0; i < ITERATIONS; i++) {
float x = curand_uniform(&rng); // Random x position in [0,1]
float y = curand_uniform(&rng); // Random y position in [0,1]
counter[threadIdx.x] += 1 - int(x * x + y * y); // Hit test
}
// The first thread in *every block* should sum the results
if (threadIdx.x == 0) {
// Reset count for this block
totals[blockIdx.x] = 0;
// Accumulate results
for (int i = 0; i < WARP_SIZE; i++) {
totals[blockIdx.x] += counter[i];
}
}
}
int main(int argc, char **argv) {
int numDev;
cudaGetDeviceCount(&numDev);
if (numDev < 1) {
cout << "CUDA device missing! Do you need to use optirun?\n";
return 1;
}
cout << "Starting simulation with " << NBLOCKS << " blocks, " << WARP_SIZE << " threads, and " << ITERATIONS << " iterations\n";
// Allocate host and device memory to store the counters
Count *hOut, *dOut;
hOut = new Count[NBLOCKS]; // Host memory
cudaMalloc(&dOut, sizeof(Count) * NBLOCKS); // Device memory
// Launch kernel
picount<<<NBLOCKS, WARP_SIZE>>>(dOut);
// Copy back memory used on device and free
cudaMemcpy(hOut, dOut, sizeof(Count) * NBLOCKS, cudaMemcpyDeviceToHost);
cudaFree(dOut);
// Compute total hits
Count total = 0;
for (int i = 0; i < NBLOCKS; i++) {
total += hOut[i];
}
Count tests = NBLOCKS * ITERATIONS * WARP_SIZE;
cout << "Approximated PI using " << tests << " random tests\n";
// Set maximum precision for decimal printing
cout.precision(DblLim::max_digits10);
cout << "PI ~= " << 4.0 * (double)total/(double)tests << endl;
return 0;
}
Il codice di per sé è abbastanza commentato, mi piacerebbe giusto sottolineare un paio di concetti e di soluzioni, per poi parlare delle prestazioni.
Prima di tutto, il numero di punti generati casualmente che andremo a creare dipende da tre parametri: WARP_SIZE, NBLOCKS e ITERATIONS.
Per capire cosa rappresentano, serve una basilare introduzione al funzionamento di CUDA e della computazione sulle GPU. Innanzitutto la funzione che andremo ad eseguire in parallelo viene chiamata kernel. Ogni kernel può essere eseguito in una grid, ogni grid contiene un certo numero di blocchi, i quali a loro volta contengono una certa quantià di thread. Che detto così vuol dire poco, ma ciascun sottoinsieme gode di proprietà interessanti.
I thread contenuti in ciascun blocco, ad esempio, condividono una porzione di memoria e possono scambiare dati tra di loro. Il warp invece è un insieme di 32 thread che vengono gestiti come unica unità dal processore. Essenzialmente quello che succederà dopo la chiamata di picount<<<NBLOCKS, WARP_SIZE>>>(dOut) è che verranno generati NBLOCKS blocchi (640 nel nostro caso, determinato dall’hardware a disposizione), ciascuno contenente WARP_SIZE thread (32 nel nostro caso). Ogni thread ripeterà la prova tante volte quante vale ITERATIONS.
Un’altro aspetto su cui vale la pena di posare la nostra attenzione è il modo con cui viene testato se il punto è nel quadrante di cerchio o meno: counter[threadIdx.x] += 1 – int(x * x + y * y). Ha attirato la mia attenzione, perché nel mio primo programma per verificare la condizione utilizzavo un if. Pensando che fosse dovuto a qualche tipo di ottimizzazione che non avevo ben capito, ho chiesto lumi ad Alessandro, e le ragioni si sono rivelate ben più profonde.
Il fatto è che quando si utilizza un if, si creano due possibili rami da seguire: quello per cui la condizione è verificata e quello per cui la condizione non è verificata. E questo è il momento in cui ho capito che la computazione parallela richiede un modo di ragionare totalmente diverso da quella tradizionale. Nella GPU infatti, quando si vengono a creare due rami, vengono prima eseguite le operazioni per il gruppo appartenente a quel ramo, mentre l’altro viene messo in pausa, e poi viene messo in pausa il primo e vengono eseguite le operazioni del secondo. Questo evidentemente comporta dei rallentamenti, avremo una parte della nostra potenza di calcolo ferma per gran parte del tempo.
Molto meglio quindi una condizione matematica, che viene eseguita allo stesso modo su tutti i thread, senza rallentare l’esecuzione.
Nel mio caso ho voluto portare ITERATIONS da 1’000’000 a 100’000’000, per un totale di 2’048’000’000’000 punti. Si legge due bilioni e quarantotto miliardi di punti. L’ultima volta, col multithreading eravamo riusciti a generare dieci miliardi di prove in circa tre minuti. Di seguito lo screenshot di questa prova:
Cinquantadue secondi e mezzo per fare quello che prima avremmo fatto in dieci ore. E il codice probabilmente potrebbe ancora essere ottimizzato. E siamo molto vicini alla sesta cifra decimale (3,141592). Spero di avervi dato un’idea della potenza di questa tecnica.
In conclusione voglio ringraziare Carmine e Alessandro per aver aver condiviso con me, e quindi permesso di condividere con chi legge, tutte queste informazioni.
E a questo punto la domanda che spero sorga spontanea è: “Ma Valerio, e tu cosa hai fatto in tutto questo? Tutto il codice che ci hai mostrato è stato scritto da altri, cosa c’è di tuo? ” Beh, in tutto questo io sono un po’ il narratore. Un narratore interno, come piaceva classificarlo ai miei insegnanti di italiano. Diciamo che io ho raccolto la storia, ne ho messo i pezzi insieme e l’ho raccontata.
Alla prossima!
This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.