Calcul Stiintific de Inalta Performanta pe Procesoare Grafice

Cuprins

Motivarea alegerii temei de licență

Tema lucrării mele de licență se intitulează „Calcul științific de înaltă performanță pe procesoare grafice” și reprezintă o curiozitate intelectuală și o nouă provocare în materie de dezvoltarea programelor.

Motivul alegerii acestei teme rezultă din curiozitatea de documentare cu privire la procesoarele grafice. Fiind o temă de actualitate, așadar un lucru nou pentru mine, am rămas surprins de folosirea unităților de procesare video (GPU) în rezolvarea problemelor complexe. Observând că în ultima perioadă de timp unitățile de procesare video sunt folosite pentru diverse calcule științifice complexe, exemplu fiind problema celor N corpuri (în tratarea clasică dar și cuantică), am rămas plăcut surprins să mă documentez cu referire la tehnologia CUDA a celor de la Nvidia.

Sunt surprins de platforma CUDA, deoarece a arătat potențial încă de la introducerea sa, ceea ce-i oferă un mare plus deoarce în momentul de față nu există o concurență acerbă de piață în materie de GPU Compute .

Procesoarele grafice au avut în ultimii ani o evoluție interesantă, apropiindu-se de statutul de co-procesor al unității centrale de prelucrare. În momentul de față plăcile video, pe lângă funcțiile specializate de prelucrare și redare a imaginilor oferă avantaje pentru realizarea calculelor matematice complexe.

În ziua de astăzi, la nivel de „Stiinta Calculatoarelor” se vorbește din ce în ce mai mult de calcul paralel, ceea ce este defapt o execuție în paralel pe mai multe procesoare a acelorași instrucțiuni, sau și a unor instrucțiuni diferite, menite în scopul rapidității rezolvării a unei probleme adaptate. Unitățile de procesare centrală și grafică aveau modul de funcționare separat, dar în recenta perioadă de timp s-a reușit integrarea CPU/GPU, ca rezultat fiind mărirea puterii de calcul cu peste 20%. Performanța a fost obținută prin conectarea celor două procesoare în așa fel încât fiecare profită de avantajele celuilalt : GPU executând calcule, iar CPU preincarcand date utilizate apoi de GPU din memorie. Prelucrarea de tip secvențial a informației a fost înlocuită cu algoritmi proiectați pentru rularea în paralel a mai multor fire de execuție.

Un alt lucru care m-a determinat să studiez această arhitectură dedicată procesării paralele este faptul că este disponibil către utilizatori suportul software pentru implementarea algoritmilor dar și suportul hardware, lângă acestea fiind atașate un set de documentații și exemple.

CAPITOLUL 1. PREZENTAREA PLATFORMEI CUDA

1.1.Introducere CUDA

CUDA este o arhitectură software și hardware pentru calcul paralel al datelor ce oferă capacități pentru procesarea cu ajutorul procesoarelor grafice denumite și GPU (graphical processing unit).

Unele dintre cele mai complexe aplicații de date existente sunt realizate combinând tehnologia de procesare a procesoarelor grafice și mediul software de lucru CUDA.

Știind că API ( Application Programming Interface ) reprezintă o interfață pentru programele de aplicații, care stabilește în amănunt modul în care programele pot accesa serviciile sistemului de operare sub care acesta rulează, CUDA oferă o serie de avantaje față de interfețele tradiționale de prelucrare a datelor cu ajutorul procesoarelor grafice. Unele dintre marile avantaje a folosirii CUDA sunt citirile nesecvențiale, memoria partajată, descărcări rapide spre procesorul video dar și un suport complet pentru operații cu numere întregi.

Procesorul grafic este considerat ca un dispozitiv de calcul având rolul de co-procesor pentru unitatea centrală de procesare, denumită și CPU (central processing unit).

CUDA este o arhitectură de calcul paralel dezvoltat de NVIDIA și este motorul de calcul din procesoarele grafice NVIDIA. CUDA vine cu un mediu software care permite dezvoltatorilor să utilizeze C ca un limbaj de programare de nivel înalt. Printre interfețele de programare a aplicațiilor pentru procesorul grafic folosind CUDA se număra C, C++, Fortran, Java Python Wrappers, DirectCompute și multe altele, așadar CUDA este conceput pentru a sprijini diferite limbaje și interfețe de programare a aplicațiilor.

1.2.Arhitectura CUDA

Arhitectura CUDA constă în foarte multe componente de bază, printre ele fiind : motoarele de procesare paralelă din interiorul procesoarelor grafice NVIDIA, suportul pentru inițializare, configurare la nivelul nucleului sistemului de operare, driverul video la nivel de utilizator și un set de funcții și instrucțiuni pentru nucleele de prelucrare paralelă.

Spre deosebire de generațiile anterioare, arhitectura CUDA a inclus o nouă oportunitate, permițându-i fiecărei unități aritmetice și logice (ALU), de pe chip să fie adunate printr-un program care intenționează să efectueze calcule de uz general. În plus, este permis unităților de execuție pe procesorul grafic scrierea și accesul la memorie, precum și accesul la memoria partajată. Toate aceste caracteristici ale arhitecturii CUDA s-au adăugat cu scopul de a creea un procesor grafic care să exceleze la calcul și la sarcinile grafice.

Legea lui Moore descrie o tendință pe termen lung în istoria mașinilor de calcul: numărul de tranzistori care pot fi plasați pe un circuit integrat se dublează aproximativ la fiecare doi ani.

O nouă provocare este dezvoltarea de aplicații software care scalează transparent paralelismul acestuia să impulsioneze creșterea numărului de nuclee pe procesoare.

Modelul de programare paralelă CUDA este conceput pentru a depăși această provocare menținând în același timp o nouă metodă de învățare pentru programatorii familiarizați cu limbaje de programe standard, cum ar fi C.

1.2.1. Tipuri de memorie CUDA

Un dispozitiv CUDA are un număr de diferite componente ale memoriei, care sunt disponibile pentru programatori.

Aceste tipuri de memorie sunt : memoria comună, memoria locală, memoria globală și memoria constant. Figura de mai jos ilustrează modul firelor de execuție în dispozitivul CUDA.

Figura 1.1. Modul firelor de execuție.

Folosim săgeată simplă, respectiv săgeată dublă, pentru a indica capacitatea de citire respectiv capacitatea de scriere. O săgeată orientată spre o componeta de memorie indica capacitatea de scriere; o săgeată orientată dinspre o componentă de memore indică capacitatea de citire. De exemplu, memoria globală și memoria constantă pot fi citite cât si scrise.

Pe un dispozitiv CUDA, mai multe nuclee pot fi invocate. Fiecare nucleu este o rețea independentă construită din unul sau mai multe blocuri de execuție. Fiecare bloc de execuție are propria memorie per bloc, care este împărțită între firele de execuție din acel bloc. Toate subiectele pot accesa scrierea sau citirea, diferite părți ale memoriei de pe dizpozitiv fiind prezentate mai jos.

În general, vom folosi gazda pentru transferul de date către și de la memoria comună și transfer către și de la memoria permanentă.

Odată ce datele sunt în memoria aparatului, firele de execuție pot citi și scrie diferite părți ale memoriei astfel:

Tabelul 1.1. Citirea și scrierea – Memorie

Tipurile de memorie prezentate în “Tabelul 1.1” prezintă o serie de avantaje, fiind descrise in prezentarea ce urmează:

Regiștri:

cea mai rapidă formă de memorie de pe multi-procesor;

este accesibilă de către firele de execuție;

durata de viață este aceeași ca și a firului de execuție;

Memoria comună :

poate fi la fel de rapidă ca regiștrii atunci când nu există conflicte de acces între firele de execuție;

accesibilă tuturor firelor de execuție dintr-un bloc;

durata de viață este aceeași ca și a blocului firelor de execuție;

Memoria globală:

poate fi de 150 de ori mai lentă decât regiștrii sau memoria comună, mai ales pentru citiri sau scrieri divergente;

accesibilă atât de către “host” cât și de “device”;

durata de viață este acceași ca și a aplicației;

Memoria locală:

este o posibilă capcană, ea aflându-se defapt în memoria globală;

în momentul în care compilatorul nu mai are loc în registri va stoca datele în memoria locală în loc să o facă în regiștri;

este accesibilă doar de către firele de execuție;

are durata de viață ca și a firelor de execuție;

Memoria constantă:

se află în memoria globală;

citește numai contextul corespunzător “device-ului”;

ascunsă;

Memoria dedicată texturilor:

se află în memoria globală;

mod special de acces la date;

1.2.2. Mărimea și lățimea benzii

Diferitele componente ale memoriei au dimensiuni și lățime de bandă diferite. Memoria globală are cea mai mică lățime de bandă, dar dimensiunea memoriei este mare.

Pentru fiecare bloc de execuție, accesul la memoria partajată este mai rapid decât la memoria locală sau la memoria constantă dar mai lentă decât registrele fiecărui fir de execuție.

Fiecare bloc de execuție are maxim 48KB memorie partajată. Registrele per fir de execuție pot stoca doar o cantitate mică de date, dar sunt cele mai rapide. Per fir memoria locală este mai lentă decât regiștrii și este utilizată pentru păstrarea datelor ce nu încap în regiștrii.

1.2.3. Variabile Automate

O variabilă automată este accesibilă numai de către firele de execuție din memoria locală. Totuși este posibil ca compilatorul să stocheze șiruri de numere în regiștrii, dacă toate accesările acestui șir se fac cu valori constante ale indicelui.

O matrice automată poate fi declarată folosind noțiunea de mărimea și lățimea benzii în felul următor : pe blocul de execuție, memoria partajată este mai rapidă decât memoria locală și memoria constantă dar mai lentă decât registrele de pe firul de execuție. O variabilă automată este accesibilă numai de către firele de execuție din memoria locală. Această varibilă automată este declarată fără calificare (fără un cuvânt cheie specific). Variabila astfel definită va fi memorată în registrul aparținând de firul de execuție respectiv și va fi accesibilă numai acelui fir.

Pentru a arăta că variabilă este partajata într-un bloc se folosește calificativul ”__shared__” de mai jos:

int autovar;

int autoarr[];

__shared__int shvar;

Calificativul ”__device__” este folosit pentru declararea locației variabilei în memorie la nivel local, iar calificativul ”__constant__” declară o variabilă constantă care se afla în memoria constantă. O variabilă constantă trebuie declarată în afara oricărei funcții, inclusiv în afara programului principal main():

__device__ int dvvar;

__constant__ int cnvar;

Referirea prin adresă(folosind pointeri) se face atribuind unei variabile declarate în comun sau în memoria globală o variabilă pointer.

float *ptr= &GlobalVar;

1.2.4. Administrarea memoriei

Administrarea memoriei se referă la modul de alocare a spațiului de memorie și transferul de date între gazdă și dispozitiv. Administrarea memoriei pe un dispozitiv CUDA este similar cu modul în care acesta se face în programarea pe unitatea centrală de programare. Este nevoie de alocarea spațiului de memorie pe gazdă, transfer date la dispozitiv folosind API, recuperarea datelor și în cele din urmă eliberarea memoriei alocate. Toate aceste sarcini sunt efectuate de dispozitivul gazdă (”host”, denumire standard pentru CPU în sintaxa CUDA).

Declarația de alocare a memoriei este dată de sintaxa:

// Declararea pointerilor
int * d_t1;
int * h_t1;

// Alocare memorie gazda
h_t1 = (int *) malloc (dimensiune);
// Alocare memorie dispozitiv
cudaMalloc ((void **) și d_t1, dimensiune);

În realitate d_t1 și h_t1 reprezintă matricea de memorie gazda.

Transferul memoriei se face odată ce spațiul de memorie este alocat, apoi putem transfera datele la procesorul grafic. Funcția cudaMemcpy invocă nucleul și copie datele de pe primul dispozitiv (”host”) la procesorul grafic (”device”). Secvența de cod corespunzătoare este prezentată mai jos :

cudaMemcpy(d_t1, h_t1, size, cudaMemcpyHostToDevice);

//Invocarea nucleului

thread_multi<<<4,N/4>>>(d_t1);

cudaMemcpy(h_t1, d_t1, size, cudaMemcpyDeviceToHost);

1.3.Calculul paralel

Calculul paralel reprezintă execuția în paralel pe mai multe procesoare a acelorași instrucțiuni, sau a unor instrucțiuni diferite, cu scopul rezolvării mai rapide a unei probleme, de obicei special adaptată sau subdivizată.

Un sistem de calcul paralel este un calculator cu mai multe procesoare care lucrează în paralel. Există multe tipuri de sisteme de calcul paralel, deosebirile dintre ele fiind în tipul de interconectare, între procesoarele componente și între procesoare și memorie.

Algoritmii paraleli se construiesc de obicei reproiectând algoritmii seriali, în direcția utilizării resurselor specifice sistemului de calcul paralel. Nu toți algoritmii seriali pot fi paralelizați. Programele cu execuție paralelă sunt de cele mai multe ori legate strâns de hardware, scopul lor fiind acela de a obține într-un mod cât mai rapid performanță.

Taxonomia lul Flynn explică clasificarea sistemelor de calcul paralel și scalar după caracteristicile instrucțiunilor și datelor.

Cea mai simplă clasificare fiind :

SIMD – o singură instrucțiune și date multiple;

MIMD – mai multe instrucțiuni și date multiple;

O altă clasificare a sistemelor de calcul paralel este bazată pe arhitectura memoriei.

Sistemele de calcul paralel cu memorie partajată dispun de procesoare multiple care pot accesa toată memoria disponibilă ca un spațiu global ce pot fi împărțite în două mari clase, în funcție de tipul de acces la memorie:

UMA – acces uniform al memoriei, în care timpii de acces la toate părțile memoriei sunt egali;

NUMA – acces neuniform al memoriei, în care timpii de acces la memorie nu sunt egali;

Sistemele de calcul paralel pot fi de asemenea clasificate și după numărul de procesoare din componența lor. Sistemele cu mii de asemenea procesoare sunt cunoscute drept “massively

parallel systems”. Restul sistemelor de calcul paralel sunt numite ”la scară mare” sau ”la scară mică”.Acesta depinde si de viteza procesorului. De exemplu un sistem de calcul paralel bazat pe un calculator personal va fi considerat unul la scară mică.

Sistemele de calcul paralel mai pot fi împărțite în sisteme multiprocesor simetrice și asimetrice, după cum toate procesoarele sunt de acelasi fel sau nu. De exemplu dacă doar un procesor poate rula codul sistemului de operare și celelalte nu au acest privilegiu, atunci sistemul este considerat ca fiind un sistem asimetric.

1.3.1. Programarea in calculul paralel

Programarea în calcul paralel cuprinde subdivizarea problemei de rezolvat, conceperea și implementare programelor corespunzătoare. Această programare se axează pe partiționarea întregii probleme de rezolvat în sarcini separate, alocarea sarcinilor procesoarelor disponibile și ulterior, sincronizarea lor pentru a obține rezultatele dorite. Acest tip de programare se poate aplica în general doar problemelor care sunt concepute în așa fel în cât pot fi paralelizabile.

În programarea paralelă există două tipuri de abordare a problemei, una dintre ele fiind paralelism implicit iar cea de-a doua fiind paralelismul explicit. Diferență este că în prima situație compilatorul descompune problema și aloca sarcini ficarui procesor în mod automat iar în calculul paralel explicit programatorul trebuie să precizeze explicit în program modul cum va fi partiționată problema de rezolvat.

1.3.2. Paralelizare la nivel de instrucțiune

Paralelizarea la nivel de instrucțiune este denumirea dată unei suite de metode de proiectare pentru majoritatea familiilor de procesoare și compilatoare în scopul măririi vitezei de execuție, mărire a vitezei de execuție ce rezultă din rularea în paralel a mai multor operații.

În spatele paralelizării la nivel de instrucțiune se afla mărirea vitezei de execuție a instrucțiunilor ce în totalitate alcătuiesc un program. Viteza de execuție se măsoară în număr de cicluri per instrucțiune.

1.3.3. Evaluarea performanțelor

Prin paralelizarea unui program secvențial se urmărește în primul rând obținerea unui timp de execuție cât mai mic ce poate fi comparat cu timpul secvențial de execuție.

Accelerarea paralelă, se calculează ca raportul dintre timpul secvențial de execuție și timpul de execuție paralelă .

(1.1)

Valoarea maximă a accelerării paralele este egală cu numărul de procesoare din sistem. O astfel de valoare poate fi atinsă într-un sistem ideal în care nu există costuri de comunicare iar procesoarele sunt echilibrat încărcate.

(1.2)

1.3.4. Legea lui Amdahl

Legea lui Amdahl definește limitările fundamentale ale programării paralele, propunând o evaluare detaliată a performanțelor mașinilor de calcul.

În legea prezentată se definește un factor de accelerare, prin care se poate stabili de câte ori mașina îmbunătățită va rula mai repede. Acesta se calculeza ca un raport dintre timpul de execuție a unei comenzi fără îmbunătățire și timpul de execuție folosind o îmbunătățire când este posibil.

(1.3)

Din timpul total de rulare, orice îmbunătățire are o fracțiune de timp de utilizare. Pe baza formulei de mai sus se poate calcula un factor de accelerare numai pentru anumite îmbunătățiri:

(1.4)

unde este fracțiunea din timpul de calcul a mașinii inițiale care poate fi convertită pentru a putea profita de îmbunătățirea făcută, iar este îmbunătățirea, ce este calculate că raportul dintre noul timp de execuție și timpul vechi de execuție, aceste două fiind la nivel de comandă.

(1.5)

În această lege se presupune că mărimea problemei care se rulează rămâne acceași când este paralelizată.

Un exemplu concret, fiind faptul că dacă pentru o problemă dată, o implementare paralelă a algoritmului de rezolvare poate rula 20% din timpul operațiilor, cealaltă parte de 80% fiind neparalelizata, legea lui Amdahl spune că accelerarea maximă a versiunii paralelizate este de 1/(1-0.20) = 1,25 , adică versiunea paralelizată este de 1,25 ori mai rapidă decât versiunea ne-paralelizată.

1.3.5. Principiul de alocare a lui Brent

Teorema lui Brent spune că dacă un calcul paralel are k faze 1, 2, . . ., k, care necesită unități de timp, și care folosesc succesiv un număr de procesoare, atunci acest calcul poate fi efectuat folosind p procesoare într-un timp :

(1.6)

iar numărul p al procesoarelor poate fi ales arbitrar.

Toate fazele la un loc necesită asadar urmatorul timp de executie :

(1.7)

1.3.6. Exemplu de paralelizare

În codul de mai jos este explicat exemplul de paralelizare pentru problema de aflare a elementului cu valoare maximă dintr-un vector.

Secvential:

Input : A,n, vectorul A de numere intregi, de lungime

Output: m=max{A[i] | i=1. . . }

Function Max(A,n)[]

For i=1…n/2 do

B[i] max(A[2i-1],A[2i]);

If n == 2 then

Return B[1]

Else

Return max( B,n/2)

Paralel:

Input : A,n, vectorul A de numere intregi, de lungime

Output: m=max{A[i] | i=1. . . }

Function Max(A,n)[]

For i=1…n/2 do

: B[i] max(A[2i-1],A[2i]); # executie paralela

If n == 2 then

Return : B[1]

Else

Return max( B,n/2)[ ]

CAPITOLUL 2. PERFORMANȚELE CPU ȘI GPU

2.1.Prezentare CPU și GPU

Unitatea centrală de prelucrare sau unitatea centrală de procesare, denumită și CPU este elementul hardware dintr-un sistem informatic care execută instrucțiunile unui program de calculator realizând operațiuni aritmetice, logice, precum și operațiunile de intrare și ieșire ale sistemului.

Procesorul grafic denumit și GPU este un circuit electronic specializat conceput pentru a manipula rapid și a modifica memoria, aîn scopul accelerării creării de imagini într-un tampon de cadre destinant ecranului de ieșire.

Un mod simplu de a înțelege diferența dintre un procesor obișnuit și un procesor grafic este de a compara modul în care acestea prelucrează sarcini. Un procesor este format din câteva nuclee optimizate pentru procesarea secvenței de tip serial, în timp ce procesorul grafic are o arhitectură paralelă formată din mii de nuclee mici, mai eficiente concepute pentru manipularea a mai multor sarcini simultan.

Procesorul grafic se folosește la accelerarea calculului, modalitate posibilă folosind unitatea de procesare grafică împreună cu un procesor. În prezent accelerarea calculului se folosește la aplicații științifice, de analiza, de inginerie, de consum și de întreprindere, oferand performante fără precedent.

Figura 2.1. Nuclee CPU – GPU

Procesoarele grafice au fost folosite inițial doar pentru randare grafică. Deoarece tehnologia a evoluat, numărul mare de nuclee în GPU față de procesoare a fost exploatat de dezvoltarea capacității de calcul, astfel încât acestea să poată procesa simulan. În timp ce GPU poate avea sute sau chiar mii de procesoare, fiecare dintre ele rulează mai lent decât un nucleu CPU și au chiar mai puține caracteristici. Caracteristicile care lipsesc de la procesoarele grafice include memoria virtuală care este necesară pentru punerea în aplicare a unui sistem de operare modern.

Cu alte cuvinte, procesoarele și GPU au arhitecturi semnificativ diferite care le fac mai potrivite pentru diferite sarcini. Un GPU poate prelucra cantități mari de date în mai multe fluxuri, efectuare operațiunilor fiind relativ simplă. Un procesor este mult mai rapid în termeni de instrucțiuni pe secundă și poate efectua mai ușor operații complexe pe un singur sau câteva fluxuri de date, dar nu se poate ocupa în mod eficient de mai multe fluxuri simultan.

Procesoarele grafice folosesc o arhitectură fundamental diferită. Programarea GPU, fiind o tehnică diferită, include noi limbaje de programare, modificări de limbi existente, și noi paradigme de programare care sunt mai bine adaptate la exprimarea unui calcul că o operațiune în paralel să fie executată de mai multe procesoare.

În ziua de astăzi, GPU-urile sunt foarte eficiente la manipularea grafică pe calculator și de procesare a imaginii, iar structura lor paralelă le face mai eficiente decât procesoarele de uz general pentru algoritmi unde prelucrarea blocurilor mari de date se face în paralel. Într-un calculator personal, un GPU poate fi prezent pe placa de bază sau, intern, în componența anumitor procesoare.

2.2. Aplicații în CUDA

CUDA C oferă o cale simplă pentru utilizatorii familiarizați cu limbajul de programare C pentru a scrie cu ușurință programe pentru executarea de către dispozitiv. Extensiile de bază permit programatorilor să definească un nucleu că o funcție C.

Aplicațiile ce folosesc tehnologia CUDA sunt încă puține, dar numărul lor crescând în fiecare zi, clusterele de procesoarele fiind înlocuite din ce în ce mai mult de plăci video ce au un cost mai mic și sunt mult mai eficiente și rapide. Deschiderea imensei puteri de procesare a acceleratoarelor grafice pentru orice fel de aplicație e un pas tehnologic foarte important, aducând un salt major de viteză.

Aplicațiile CUDA iau amploare într-un ritm foarte bun gândindu-ne că deja există aplicații destinate simulărilor meteorologice, previziunilor financiare ce lucrează pe placi video cu o eficiență de remarcat.

În aplicațiile de mai jos, am făcut o analiză în scopul determinării timpului de execuție a rulării programului arătând eficienta efectuării calculelor exacte utilizând bibliotecile inițiale, în prima instanță pe GPU comparativ cu cea de pe CPU .

2.2.1.Înmulțirea matricilor

O matrice este o structură de numere ce ne-o imaginăm dreptunghiulară și nimic mai mult. Dar, în ciuda simplității ei, este una dintre obiectele matematice cele mai utile și fundamentale în calculul științific.

Înmulțirea matricilor este o operație fundamentală pentru calculul științific. Mai mult decât atât, modelele algoritmice de înmulțire a matricilor sunt reprezentantive. Fiind una dintre cele mai importante programe de învățare a programării paralele în continuare voi demonstra o eficientă a calculului folosind CPU și GPU.

Luăm în considerare cazul cel mai simplu de înmulțire, cu două matrici inițiale A și B, iar rezultatul îl reținem într-o matrice C.

Implementarea pe CPU este reprezentată de codul :

void main(){ 

define A, B, C

for i = 0 to M do

for j = 0 to N do

/* calculul elementului C(i,j) */

for k = 0 to K do

C(i,j) <= C(i,j) + A(i,k) * B(k,j)

end

end

end

}

Pentru a simplifica explicația, considerăm matricile pătratice ( M == N == K ) utilizate în ilustrație. Figura (2.2) arata amprenta de memorie pentru a calcula elementul C(3,11) reprezentat de culoarea roșie. Acest lucru poate fi privit că produsul intern al unui rând a lui A (în albastru) cu o coloană a lui B (în verde).

Figura 2.2. Calcul înmulțire matrice pe CPU

Implementarea pe GPU este reprezentată de codul :

void main(){

define A_cpu, B_cpu, C_cpu in the CPU memory

define A_gpu, B_gpu, C_gpu in the GPU memory

memcopy A_cpu to A_gpu

memcopy B_cpu to B_gpu

dim3 dimBlock(16, 16)

dim3 dimGrid(N/dimBlock.x, M/dimBlock.y)

matrixMul<<<dimGrid, dimBlock>>>(A_gpu,B_gpu,C_gpu,K)

memcopy C_gpu to C_cpu

}

 __global__ void matrixMul(A_gpu,B_gpu,C_gpu,K){

temp <= 0

i <= blockIdx.y * blockDim.y + threadIdx.y // Randul I al matricii C

j <= blockIdx.x * blockDim.x + threadIdx.x // Coloana j a matricii C

for k = 0 to K-1 do

accu <= accu + A_gpu(i,k) * B_gpu(k,j)

end

C_gpu(i,j) <= accu

}

Figura 2.3. Calcul înmulțire matrice pe GPU

Implementarea pe GPU atribuie unui fir de execuție sarcina de a calcula un element al matricei C. Fiecare fir conține un rând de matrice A și o coloană a matricei B din memoria globală, făcând produsul intern și apoi trimițând rezultatul înapoi matricei C în memoria globală.

Pentru a se observa diferența între timpul de execuție pe CPU și timpul de execuție pe GPU, în antetul programului se atașează un cod pentru cronometrarea timpului de execuție, urmând ca datele rezultate să le folosim pentru reprezentare.

În problema „Înmulțirea Mătricilor” se observă că timpii de execuție diferă foarte mult, plecând de la o matrice pătratica mică și ajungând la o matrice de ordin mare.

În tabelul de mai jos sunt reprezentate datele acumulate folosind cronometru, apoi reprezentarea acestora pe grafic.

Tabelul 2.1. Date acumulate pentru problema Înmulțirea Matricilor

Figura 2.4. Reprezentarea timpului de execuție a programului folosind CPU, respectiv GPU

2.2.2. Permutari CUDA C

În aplicația “Permutări CUDA C” programul este scris așa încât să lase utilizatorul să introducă de la tastatură numărul elementelor, acesta realizând automat toate permutările posibile aferente numărului de elemente.

Ca și aplicația “Înmulțirea Mătricilor” folosim un cod pentru cronometrarea timpului de execuție.

Tabelul 2.2. Date acumulate pentru problema Permutări CUDA C

Figura 2.5. Reprezentarea timpului de execuție a programului folosind CPU, respectiv GPU

CAPITOLUL 3. SIMULARE RAPIDĂ A PROBLEMEI CELOR N CORPURI CU AJUTORUL CUDA

3.1.Introducere

O simulare N corpuri aproximează numeric evoluția unui sistem de corpuri în care fiecare corp interacționează în mod continu cu alt corp.

Un exemplu familiar este o simulare de astrofizică în care fiecare corp reprezintă o galaxie, iar corpurile se atrag reciproc prin forța gravitațională.

Figura 3.1. Fereastră de la o randare interactivă 3D a unui sistem de 10240 corpuri, simulate prin aplicație

Codul sursă corespunzator programului analizat este prezentat in “Anexa 1”.

Turbulența fluxului de curgere a fluidelor și calculul de iluminare la nivel mondial în grafica pe calculator sunt alte exemple de probleme care folosesc simularea a N corpuri.

Simulări prin considerarea interacțiunile intre perechile este o tehnică ce folosește o modalitate accesibilă ce evaluează toate interacțiunile în perechi dintre corpuri. Este o metodă relative simplă dar care nu este utilizată în general în simularea sistemelor mari, din cauza complexității computaționale.

Interacțiunile între perechi sunt de obicei folosite doar pentru determinarea forțelor de interacțiune de prim rang. Metoda poate fi apoi combinată cu o metodă mai rapidă bazată pe o aproximare a câmpurilor de forțe cu rază lungă de acțiune, care este valabilă doar între părțile sistemului ce sunt bine separate.

Parcurgerea tuturor algoritmilor menționați mai sus necesită un timp substanțial pentru calcul și este prin urmare un obiectiv interesant pentru accelerare. Îmbunătățirea performanței unui algoritm ce consideră toate perechile va îmbunătăți de asemenea și performanța componentei de domeniu depărtat, de câmp mediu.

În acest capitol, ne vom concentra pe nucleul de calcul, pe perechile de calcul și pe punerea în aplicare a algoritmului în conformitate cu modelul de programare CUDA NVIDIA.

3.2.Perechile de simulare “N corpuri”

Folosim potențialul gravitațional pentru a ilustra forma de bază de calcul a tuturor perechilor de simulare N corpuri. În următorul calcul, vom folosi caractere îngroșate pentru a ne referi la mărimi vectoriale.. Având în vedere N corpuri cu poziții inițiale xi si viteze vi pentru , vectorul forța f ij care acționează asupra corpului i cauzat de atracția gravitațională a corpului j este dată in formula :

(3.1)

unde si sunt masele corpurilor i si j, respectiv este vectorul de poziție de la corpul i la corpul j iar G este constanta gravitațională.

Factorul din stânga, mărimea forței, este proporțional cu produsul maselor și scade cu pătratul distanței dintre corpurile i și j. Factorul din dreapta dă direcția forței, un vector unitate în direcția de la corpul i la corpul j (pentru că gravitația este o forță de atracție).

Forța totală F i ce acționează asupra corpului i, ca urmare a interacțiunilor sale cu celelalte N-1 corpuri, este obținută prin insumarea tuturori interacțiunilor:

(3.2)

Fiind corpuri apropiate unul de altul, forța între ele poate crește fără limită, ceea ce duce la o situatie nedorită pentru integrarea numerică. În simulările astrofizice, coliziunile dintre corpuri sunt în general excluse. Acest lucru este rezonabil dacă corpurile reprezintă galaxii sausisteme planetare, care pot trece una prin cealaltă fără coliziuni importante ale corpurilor ce le compun. Pentru a simula acest aspect, se adaugă un factor de dedurizare 2 > 0 , iar numitorul este rescris dupa cum urmează:

(3.3)

Condiția nu mai este necesara in sumă , pentru că fii = 0 iar atunci 2 > 0 . De fapt factorul de dedurizare (denumit si factorul de înmuiere) limitează mărimea forței între corpuri, de dorit pentru integrarea numerică a stării sistemului.

Pentru a integra în timp mișcarea sistemului, avem nevoie de accelerația ai = Fi /mi, pentru găsirea pozitiei si vitezei fiecarui corp i.

(3.4)

Alegerea metodei de integrare în problemele N corpuri depinde de obicei de natura sistemului studiat.

3.3.Implementarea folosind CUDA a algoritmului de calcul a interacțiunilor dintre toate perechile unui sistem de N corpuri

Ne putem gândi la algoritmul tuturori perechilor, ca la calcularea fiecărei intrări fij intr-o retea NxN a tuturor forțelor de perechi. Atunci forța totală Fi (sau acceleratia ai) a corpului i se obține din suma tuturor intrarilor in randul i. Fiecare intrare poate fi calculată independent.

Cu toate acestea, aceasta abordare necesita memorie si ar fi limitata in mod substantial de latimea de banda. In schimb am serializat unele calcule pentru a realiza reutilizarea datelor necesare pentru a atinge performanta maxima a unitatilor aritmetice si a reduce latimea de banda de memorie necesara.

Prin urmare, vom introduce noțiunea de “fragment de calcul”, o regiunea pătrată a grilei de perechi de forțe formate din p rânduri și p coloane. Doar 2p descrieri ale corpului sunt necesare pentru a evalua toate p2 interacțiuni din fragment (care pot fi refolosite ulterior). Aceste descrieri ale corpului pot fi stocate în memoria partajată sau în regiștrii. Efectul total al interacțiunilor este capturat ca o actualizare a vectorilor de accelerație p.

Pentru a realiza reutilizarea optimă a datelor am aranjat calculul, astfel încât interacțiunile din fiecare rând sunt evaluate în ordine secvențială, actualizarea vectorului de accelerație, în timp ce rândurile distincte sunt evaluate în paralel.

În figura de mai jos, diagrama din stânga arată strategia de evaluare iar în diagrama din partea dreaptă este prezentată intrările sau ieșirile calculului din “fragment”.

Figura 3.2. Figura schematica a unei placi de calcul

3.3.1. Forța de calcul corp-corp

Interacțiunea dintre o pereche de corpuri, descrisă în secțiunea anterioară este implementată ca un calcul complet în serie. Codul de mai jos calculează forță corpului i la interacțiunea cu corpul j în urmă accelelrarii a corpului i că urmare a acestei interacțiuni.

Există 20 de operații cu virgulă mobilă în acest cod, de numărare, adunare, multiplicări, apel la sqrt (radical) etc.

   __device__ float3  

bodyBodyInteraction(float4 bi, float4 bj, float3 ai)  

{  

  float3 r;  

  // r_ij  [3 FLOPS]  

  r.x = bj.x – bi.x;  

  r.y = bj.y – bi.y;  

  r.z = bj.z – bi.z;  

  // distSqr = dot(r_ij, r_ij) + EPS^2  [6 FLOPS]  

   float distSqr = r.x * r.x + r.y * r.y + r.z * r.z + EPS2;  

  // invDistCube =1/distSqr^(3/2)  [4 FLOPS (2 mul, 1 sqrt, 1 inv)]  

   float distSixth = distSqr * distSqr * distSqr;  

  float invDistCube = 1.0f/sqrtf(distSixth);  

  // s = m_j * invDistCube [1 FLOP]  

   float s = bj.w * invDistCube;  

  // a_i =  a_i + s * r_ij [6 FLOPS]  

  ai.x += r.x * s;  

  ai.y += r.y * s;  

  ai.z += r.z * s;  

  return ai;  }  

Folosim tipul de date “float4 CUDA” pentru descrierea corpului și a accelerațiilor stocate în memoria GPU. Vom stoca masa fiecărui corp și poziția corpului. Utilizarea float4 permite accesul la memoria unitate.

3.3.2. Fragmente de calcul

O piesă este evaluată de p fire de execuție ce efectuează acceași secvență de operații dar pe date diferite. Fiecare fir actualizează accelerația unui singur corp, ca urmare a interacțiunilor sale toate celelalte corpuri. Am încărcat descrierea corpului p din memoria GPU în memoria partajată oferită fiecărui bloc de fire de execuție în modemul CUDA. Fiecare fir evaluează doar interacțiuni cu indici succesivi.

Într-un bloc de fire de execuție fiecare spațiu este dimensionat pentru a echilibra paralelismul cu reutilizarea datelor. Gradul de paralelism este dat de numărul de rânduri și trebuie să fie suficient de mare pentru a utiliza cât mai bine resursele de calcul. Gradul de utilizare a datelor crește odată cu numărul de coloane, iar acest parametru reglementează, de asemenea și dimensiunea transferului datelor caracteristice unor corpuri din memoria calculatorului în memoria partajată. În final, dimensiunea fragmentelor determină și spațiul regiștrilor.

Pentru această implementare, am folosit fragmente pătratice de dimensiuni “p pe p”. Înainte de a executa un fragment de calcul, fiecare fir preia informația în memoria partajată, după care firele sunt sincronizate. Prin urmare, fiecare placă de calcul începe cu corpurile succesive “p” în memoria partajată.

Figura 3.3 arată un bloc de fire de execuție care execută cod pentru mai multe date. Timpul este prezentat pe direcție orizontală în timp ce paralelismul pe direcție verticală. Liniile continue delimitează plăcile de calcul, arătând unde memoria partajată este încărcată și unde o sincronizare este realizată. Într-un bloc de fire de execuție există N/p plăci. Fiecare fir calculează toate interacțiunile N pentru un corp.

Figura 3.3. Exemplificare bloc fire de executie – placi calcul

În exemplul de mai jos este prezentat evaluarea interacțiunilor într-o placă de calcul “pxp”:

__device__ float3

tile_calculation(float4 myPosition, float3 accel)

{

int i;

extern __shared__ float4[] shPosition;

for (i = 0; i < blockDim.x; i++) {

accel = bodyBodyInteraction(myPosition, shPosition[i], accel);

}

return accel;

}

Codul pentru a calculul forțelor dintr-un bloc de fire de execuție este prezentat mai jos. Acest cod este în fapt “kernel-ul CUDA” care este numit și gazdă.

Parametrii funcției “calculate_forces()” sunt pointeri la memoria globală cu pozițiile devX și accelerațiile devA ale corpurilor. Noi le atribuim indicii locali cu conversie de tip astfel încât să poată fi indexați ca matrici. În cazul fragmentelor de calcul este nevoie de două sincronizări, prima asigurând ca toate locațiile de memorie partajată sunt populate înainte de aplicării calculului forței gravitaționale, iar a două se asigură că toate firele au finalizat calculul lor de forță gravitațională înaintea avansării la fragmentul următor. Fără aceasta a două sincronizare, firele de execuție care termină partea lor în calcul ar putea suprascrie memoria partajată.

Kernel-ul CUDA executat de un bloc de fire de execuție cu p fire pentru a calcula accelerația gravitațională pentru p corpuri ca rezultat a tuturor N interacțiuni are sintaxa de mai jos :

__global__ void

calculate_forces(void *devX, void *devA)

{

extern __shared__ float4[] shPosition;

float4 *globalX = (float4 *)devX;

float4 *globalA = (float4 *)devA;

float4 myPosition;

int i, tile;

float3 acc = {0.0f, 0.0f, 0.0f};

int gtid = blockIdx.x * blockDim.x + threadIdx.x;

myPosition = globalX[gtid];

for (i = 0, tile = 0; i < N; i += p, tile++) {

int idx = tile * blockDim.x + threadIdx.x;

shPosition[threadIdx.x] = globalX[idx];

__syncthreads();

acc = tile_calculation(myPosition, acc);

__syncthreads();

}

// Salvare rezultat în memorie la nivel global pentru etapa de integrare.

float4 acc4 = {acc.x, acc.y, acc.z, 0.0f};

globalA[gtid] = acc4;

}

3.3.3. Definirea grilei blocurilor firelor de executie

Programul “kernel” de mai sus calculează accelerația a p corpuri din sistem, cauzată de interacțiunea lor cu toate N corpurile din sistem. Invocăm acest program pe o grilă a blocurilor firelor de execuție pentru a calcula accelerarea tuturor corpurilor N. Pentru că există p fire per bloc de execuție și un fir pe corp, numărul de blocuri necesare pentru a finaliza toate corpurile N este N/p, așa că am definit o grilă 1D de mărime N/p.

Rezultatul este un total de N fire de execuție care calculează N forțe pentru un total de interactiuni.

Evaluarea grilei complete a interacțiunilor poate fi vizualizată așa cum este prezentată în figură 3.4. Axa verticală prezintă paralelismul grilei 1D, N/p blocuri independente cu p fire fiecare. Axa orizontală reprezintă prelucrarea secvențială a N calcule în fiecare fir. Un bloc de fire de execuție încărca memoria partajată în p păși pentru a partaja datele.

Figura 3.4. Grila de blocuri care calculează toate forțele

3.4.Rezultate de performanță

Întorcăndu-ne la simularea noastră a sistemului de 10240 corpuri prin simplă vizualizare a capacității GeForce GT 730M (placă folosită în cazul problemelor rezolvate) observăm că aceasta este capabilă de generarea a 10240 corpuri, având loc calculul a 9410 miliarde de interacțiuni pe secundă.

Figura 3.5. Imagine de la rularea programului

În schimb, instrucțiuni complexe, cum ar fi rădăcina pătrată au nevoie de mai multe cicluri de procesor.

Atunci când se compară ratele gigaflop, pur și simplu se număra operațiunile virgulă mobilă. Astfel, prin numărarea operațiunilor în virgula mobilă ( în codul prezentat în capitolul “Fortele de calcul corp-corp”), identifica, adunări și înmulțiri, un calcul al rădăcinei pătrate și o împărțire. Împărțirea și rădăcina pătrată necesită în mod evident mai mult timp decât operațiile de adunare sau înmulțire.

Creșterea performanței este subdivizata în 3 categorii:

cu derulare in bucla;

cand dimensiunea blocului de fire de execuție variază;

pentru N corpuri mici;

Figura 3.6. Cresterea performantei folosind multiple fire de executie per corp

Figura 3.7. Fereastra de la o randare interactiva 3D a unui sistem de 10240 corpuri, simulate prin aplicatie folosind precizie dubla

3.5.Concluzie

Este dificil să ne imaginăm un algoritm al lumii reale, care este mai bine adaptat la executarea pe arhitectura CUDA decât algoritmul ce calculează interacțiunile între toate perechile dintr-un sistem de N corpuri. În acest capitol am demonstrat trei caracteristici ale algoritmului care ajută să se realizeze o astfel de eficiență:

paralelism simplu cu modele de acces la memoria secventiala;

reutilizarea datelor ce tin ocupate unitatile aritmetice;

operatiuni complexe cum ar fi inversul radacinii patrate;

Rezultatul este un algoritm care rulează de aproximativ 50 de ori mai repede decât o implementare serială

La acest nivel de performanța, simulări 3D cu un număr mare de particule pot fi rulate interactiv.

CUDA poate asigură viteze mari pentru executarea datelor în paralel, viteze și mai mari fiind posibile cu un hardware mai bun și o administrare corespunzătoare.

Știind că ”FLOPS” reprezintă o unitate de măsură a puterii de calcul a unui calculator măsurând numărul maxim de operații în virgule mobile ce sunt executate pe secundă, o performanță deosebită o au calculele matematice care ajung la 472 GFLOPS și o bandă de memorie de 80 GB/s. În cazul performanței problemei N corpuri, vizualizând Figura 3.6 se poate vedea o performanță considerabilă făcând analogia cu firele de execuție corespunzătoare programului.

ANEXA 1. Codul programului “N corpuri”

“bodysystem.h”

#ifndef __BODYSYSTEM_H__

#define __BODYSYSTEM_H__

#include <algorithm>

enum NBodyConfig

{

NBODY_CONFIG_RANDOM,

NBODY_CONFIG_SHELL,

NBODY_CONFIG_EXPAND,

NBODY_NUM_CONFIGS

};

enum BodyArray

{

BODYSYSTEM_POSITION,

BODYSYSTEM_VELOCITY,

};

template <typename T> struct vec3

{

typedef float Type;

}; // model

template <> struct vec3<float>

{

typedef float3 Type;

};

template <> struct vec3<double>

{

typedef double3 Type;

};

template <typename T> struct vec4

{

typedef float Type;

}; // model

template <> struct vec4<float>

{

typedef float4 Type;

};

template <> struct vec4<double>

{

typedef double4 Type;

};

class string;

// Clasa BodySystem abstract

template <typename T>

class BodySystem

{

public: // metode

BodySystem(int numBodies) {}

virtual ~BodySystem() {}

virtual void loadTipsyFile(const std::string &filename) = 0;

virtual void update(T deltaTime) = 0;

virtual void setSoftening(T softening) = 0;

virtual void setDamping(T damping) = 0;

virtual T *getArray(BodyArray array) = 0;

virtual void setArray(BodyArray array, const T *data) = 0;

virtual unsigned int getCurrentReadBuffer() const = 0;

virtual unsigned int getNumBodies() const = 0;

virtual void synchronizeThreads() const {};

protected: // metode

BodySystem() {} // constructor

virtual void _initialize(int numBodies) = 0;

virtual void _finalize() = 0;

};

inline float3

scalevec(float3 &vector, float scalar)

{

float3 rt = vector;

rt.x *= scalar;

rt.y *= scalar;

rt.z *= scalar;

return rt;

}

inline float

normalize(float3 &vector)

{

float dist = sqrtf(vector.x*vector.x + vector.y*vector.y + vector.z*vector.z);

if (dist > 1e-6)

{

vector.x /= dist;

vector.y /= dist;

vector.z /= dist;

}

return dist;

}

inline float

dot(float3 v0, float3 v1)

{

return v0.x*v1.x+v0.y*v1.y+v0.z*v1.z;

}

inline float3

cross(float3 v0, float3 v1)

{

float3 rt;

rt.x = v0.y*v1.z-v0.z*v1.y;

rt.y = v0.z*v1.x-v0.x*v1.z;

rt.z = v0.x*v1.y-v0.y*v1.x;

return rt;

}

// functia pentru utilitate

template <typename T>

void randomizeBodies(NBodyConfig config, T *pos, T *vel, float *color, float clusterScale,

float velocityScale, int numBodies, bool vec4vel)

{

switch (config)

{

default:

case NBODY_CONFIG_RANDOM:

{

float scale = clusterScale * std::max<float>(1.0f, numBodies / (1024.0f));

float vscale = velocityScale * scale;

int p = 0, v = 0;

int i = 0;

while (i < numBodies)

{

float3 point;

//const int scale = 16;

point.x = rand() / (float) RAND_MAX * 2 – 1;

point.y = rand() / (float) RAND_MAX * 2 – 1;

point.z = rand() / (float) RAND_MAX * 2 – 1;

float lenSqr = dot(point, point);

if (lenSqr > 1)

continue;

float3 velocity;

velocity.x = rand() / (float) RAND_MAX * 2 – 1;

velocity.y = rand() / (float) RAND_MAX * 2 – 1;

velocity.z = rand() / (float) RAND_MAX * 2 – 1;

lenSqr = dot(velocity, velocity);

if (lenSqr > 1)

continue;

pos[p++] = point.x * scale; // pozitie x

pos[p++] = point.y * scale; // pozitie y

pos[p++] = point.z * scale; // pozitie z

pos[p++] = 1.0f; // mass

vel[v++] = velocity.x * vscale; // pozitie x

vel[v++] = velocity.y * vscale; // pozitie x

vel[v++] = velocity.z * vscale; // pozitie x

if (vec4vel) vel[v++] = 1.0f; // inverse mass

i++;

}

}

break;

case NBODY_CONFIG_SHELL:

{

float scale = clusterScale;

float vscale = scale * velocityScale;

float inner = 2.5f * scale;

float outer = 4.0f * scale;

int p = 0, v=0;

int i = 0;

while (i < numBodies)//for(int i=0; i < numBodies; i++)

{

float x, y, z;

x = rand() / (float) RAND_MAX * 2 – 1;

y = rand() / (float) RAND_MAX * 2 – 1;

z = rand() / (float) RAND_MAX * 2 – 1;

float3 point = {x, y, z};

float len = normalize(point);

if (len > 1)

continue;

pos[p++] = point.x * (inner + (outer – inner) * rand() / (float) RAND_MAX);

pos[p++] = point.y * (inner + (outer – inner) * rand() / (float) RAND_MAX);

pos[p++] = point.z * (inner + (outer – inner) * rand() / (float) RAND_MAX);

pos[p++] = 1.0f;

x = 0.0f; // * (rand() / (float) RAND_MAX * 2 – 1);

y = 0.0f; // * (rand() / (float) RAND_MAX * 2 – 1);

z = 1.0f; // * (rand() / (float) RAND_MAX * 2 – 1);

float3 axis = {x, y, z};

normalize(axis);

if (1 – dot(point, axis) < 1e-6)

{

axis.x = point.y;

axis.y = point.x;

normalize(axis);

}

//if (point.y < 0) axis = scalevec(axis, -1);

float3 vv = {(float)pos[4*i], (float)pos[4*i+1], (float)pos[4*i+2]};

vv = cross(vv, axis);

vel[v++] = vv.x * vscale;

vel[v++] = vv.y * vscale;

vel[v++] = vv.z * vscale;

if (vec4vel) vel[v++] = 1.0f;

i++;

}

}

break;

case NBODY_CONFIG_EXPAND:

{

float scale = clusterScale * numBodies / (1024.f);

if (scale < 1.0f)

scale = clusterScale;

float vscale = scale * velocityScale;

int p = 0, v = 0;

for (int i=0; i < numBodies;)

{

float3 point;

point.x = rand() / (float) RAND_MAX * 2 – 1;

point.y = rand() / (float) RAND_MAX * 2 – 1;

point.z = rand() / (float) RAND_MAX * 2 – 1;

float lenSqr = dot(point, point);

if (lenSqr > 1)

continue;

pos[p++] = point.x * scale; // pos.x

pos[p++] = point.y * scale; // pos.y

pos[p++] = point.z * scale; // pos.z

pos[p++] = 1.0f; // mass

vel[v++] = point.x * vscale; // pos.x

vel[v++] = point.y * vscale; // pos.x

vel[v++] = point.z * vscale; // pos.x

if (vec4vel) vel[v++] = 1.0f; // inverse mass

i++;

}

}

break;

}

if (color)

{

int v = 0;

for (int i=0; i < numBodies; i++)

{

//const int scale = 16;

color[v++] = rand() / (float) RAND_MAX;

color[v++] = rand() / (float) RAND_MAX;

color[v++] = rand() / (float) RAND_MAX;

color[v++] = 1.0f;

}

}

}

#endif // __BODYSYSTEM_H__

“bodysistemcuda.cu”

#include <helper_cuda.h>

#include <math.h>

#include <GL/glew.h>

#if defined(__APPLE__) || defined(MACOSX)

#pragma clang diagnostic ignored "-Wdeprecated-declarations"

#include <GLUT/glut.h>

#else

#include <GL/freeglut.h>

#endif

// CUDA librarii

#include <cuda_runtime.h>

#include <cuda_gl_interop.h>

#include "bodysystem.h"

__constant__ float softeningSquared;

__constant__ double softeningSquared_fp64;

cudaError_t setSofteningSquared(float softeningSq)

{

return cudaMemcpyToSymbol(softeningSquared,

&softeningSq,

sizeof(float), 0,

cudaMemcpyHostToDevice);

}

cudaError_t setSofteningSquared(double softeningSq)

{

return cudaMemcpyToSymbol(softeningSquared_fp64,

&softeningSq,

sizeof(double), 0,

cudaMemcpyHostToDevice);

}

template<class T>

struct SharedMemory

{

__device__ inline operator T *()

{

extern __shared__ int __smem[];

return (T *)__smem;

}

__device__ inline operator const T *() const

{

extern __shared__ int __smem[];

return (T *)__smem;

}

};

template<typename T>

__device__ T rsqrt_T(T x)

{

return rsqrt(x);

}

template<>

__device__ float rsqrt_T<float>(float x)

{

return rsqrtf(x);

}

template<>

__device__ double rsqrt_T<double>(double x)

{

return rsqrt(x);

}

// Macro-uri pentru a simplifica adresarea la memorie

#define SX(i) sharedPos[i+blockDim.x*threadIdx.y]

// Acest macro este folosit doar cand ‘multithreadBodies’ este adevarat

#define SX_SUM(i,j) sharedPos[i+blockDim.x*j]

template <typename T>

__device__ T getSofteningSquared()

{

return softeningSquared;

}

template <>

__device__ double getSofteningSquared<double>()

{

return softeningSquared_fp64;

}

template <typename T>

struct DeviceData

{

T *dPos[2]; // indicii gazda

T *dVel;

cudaEvent_t event;

unsigned int offset;

unsigned int numBodies;

};

template <typename T>

__device__ typename vec3<T>::Type

bodyBodyInteraction(typename vec3<T>::Type ai,

typename vec4<T>::Type bi,

typename vec4<T>::Type bj)

{

typename vec3<T>::Type r;

// r_ij [3 FLOPS]

r.x = bj.x – bi.x;

r.y = bj.y – bi.y;

r.z = bj.z – bi.z;

// distSqr = dot(r_ij, r_ij) + EPS^2 [6 FLOPS]

T distSqr = r.x * r.x + r.y * r.y + r.z * r.z;

distSqr += getSofteningSquared<T>();

// invDistCube =1/distSqr^(3/2) [4 FLOPS (2 mul, 1 sqrt, 1 inv)]

T invDist = rsqrt_T(distSqr);

T invDistCube = invDist * invDist * invDist;

// s = m_j * invDistCube [1 FLOP]

T s = bj.w * invDistCube;

// a_i = a_i + s * r_ij [6 FLOPS]

ai.x += r.x * s;

ai.y += r.y * s;

ai.z += r.z * s;

return ai;

}

template <typename T>

__device__ typename vec3<T>::Type

computeBodyAccel(typename vec4<T>::Type bodyPos,

typename vec4<T>::Type *positions,

int numTiles)

{

typename vec4<T>::Type *sharedPos = SharedMemory<typename vec4<T>::Type>();

typename vec3<T>::Type acc = {0.0f, 0.0f, 0.0f};

for (int tile = 0; tile < numTiles; tile++)

{

sharedPos[threadIdx.x] = positions[tile * blockDim.x + threadIdx.x];

__syncthreads();

// aici se foloseste prezentarea placilor de calcul prezentate in capitolul 3

#pragma unroll 128

for (unsigned int counter = 0; counter < blockDim.x; counter++)

{

acc = bodyBodyInteraction<T>(acc, bodyPos, sharedPos[counter]);

}

__syncthreads();

}

return acc;

}

template<typename T>

__global__ void

integrateBodies(typename vec4<T>::Type *__restrict__ newPos,

typename vec4<T>::Type *__restrict__ oldPos,

typename vec4<T>::Type *vel,

unsigned int deviceOffset, unsigned int deviceNumBodies,

float deltaTime, float damping, int numTiles)

{

int index = blockIdx.x * blockDim.x + threadIdx.x;

if (index >= deviceNumBodies)

{

return;

}

typename vec4<T>::Type position = oldPos[deviceOffset + index];

typename vec3<T>::Type accel = computeBodyAccel<T>(position,

oldPos,

numTiles);

// acceleratie = forta/masa

// noua viteza = vechea viteza + acceleratia * deltaTIMP

typename vec4<T>::Type velocity = vel[deviceOffset + index];

velocity.x += accel.x * deltaTime;

velocity.y += accel.y * deltaTime;

velocity.z += accel.z * deltaTime;

velocity.x *= damping;

velocity.y *= damping;

velocity.z *= damping;

// noua pozitie = vechea pozitie + viteza * deltaTIMP

position.x += velocity.x * deltaTime;

position.y += velocity.y * deltaTime;

position.z += velocity.z * deltaTime;

// stocarea noii pozitii si viteze

newPos[deviceOffset + index] = position;

vel[deviceOffset + index] = velocity;

}

template <typename T>

void integrateNbodySystem(DeviceData<T> *deviceData,

cudaGraphicsResource **pgres,

unsigned int currentRead,

float deltaTime,

float damping,

unsigned int numBodies,

unsigned int numDevices,

int blockSize,

bool bUsePBO)

{

if (bUsePBO)

{

checkCudaErrors(cudaGraphicsResourceSetMapFlags(pgres[currentRead], cudaGraphicsMapFlagsReadOnly));

checkCudaErrors(cudaGraphicsResourceSetMapFlags(pgres[1-currentRead], cudaGraphicsMapFlagsWriteDiscard));

checkCudaErrors(cudaGraphicsMapResources(2, pgres, 0));

size_t bytes;

checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void **)&(deviceData[0].dPos[currentRead]), &bytes, pgres[currentRead]));

checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void **)&(deviceData[0].dPos[1-currentRead]), &bytes, pgres[1-currentRead]));

}

for (unsigned int dev = 0; dev != numDevices; dev++)

{

if (numDevices > 1)

{

cudaSetDevice(dev);

}

int numBlocks = (deviceData[dev].numBodies + blockSize-1) / blockSize;

int numTiles = (numBodies + blockSize – 1) / blockSize;

int sharedMemSize = blockSize * 4 * sizeof(T); // 4 floats for pos

integrateBodies<T><<< numBlocks, blockSize, sharedMemSize >>>

((typename vec4<T>::Type *)deviceData[dev].dPos[1-currentRead],

(typename vec4<T>::Type *)deviceData[dev].dPos[currentRead],

(typename vec4<T>::Type *)deviceData[dev].dVel,

deviceData[dev].offset, deviceData[dev].numBodies,

deltaTime, damping, numTiles);

if (numDevices > 1)

{

checkCudaErrors(cudaEventRecord(deviceData[dev].event));

// RULAREA PE VERSIUNI MAI VECHI DE KERNEL;

cudaStreamQuery(0);

}

// verificati daca incercarea kernelului a generat vreo eroare

getLastCudaError("Kernel execution failed");

}

if (numDevices > 1)

{

for (unsigned int dev = 0; dev < numDevices; dev++)

{

checkCudaErrors(cudaEventSynchronize(deviceData[dev].event));

}

}

if (bUsePBO)

{

checkCudaErrors(cudaGraphicsUnmapResources(2, pgres, 0));

}

}

// gener. cod

template void integrateNbodySystem<float>(DeviceData<float> *deviceData,

cudaGraphicsResource **pgres,

unsigned int currentRead,

float deltaTime,

float damping,

unsigned int numBodies,

unsigned int numDevices,

int blockSize,

bool bUsePBO);

template void integrateNbodySystem<double>(DeviceData<double> *deviceData,

cudaGraphicsResource **pgres,

unsigned int currentRead,

float deltaTime,

float damping,

unsigned int numBodies,

unsigned int numDevices,

int blockSize,

bool bUsePBO);

“nbody.cpp”

#include <GL/glew.h>

#if defined(WIN32) || defined(_WIN32) || defined(WIN64) || defined(_WIN64)

#include <GL/wglew.h>

#endif

#if defined(__APPLE__) || defined(MACOSX)

#pragma clang diagnostic ignored "-Wdeprecated-declarations"

#include <GLUT/glut.h>

#else

#include <GL/freeglut.h>

#endif

#include <paramgl.h>

#include <cstdlib>

#include <cstdio>

#include <algorithm>

#include <assert.h>

#include <math.h>

#include <cuda_runtime.h>

#include <cuda_gl_interop.h>

#include <helper_cuda.h>

#include <helper_cuda_gl.h>

#include <helper_functions.h>

#include "bodysystemcuda.h"

#include "bodysystemcpu.h"

#include "render_particles.h"

#include "cuda_runtime.h"

// vizualizare parametrii

int ox = 0, oy = 0;

int buttonState = 0;

float camera_trans[] = {0, -2, -150};

float camera_rot[] = {0, 0, 0};

float camera_trans_lag[] = {0, -2, -150};

float camera_rot_lag[] = {0, 0, 0};

const float inertia = 0.1f;

ParticleRenderer::DisplayMode displayMode =

ParticleRenderer::PARTICLE_SPRITES_COLOR;

bool benchmark = false;

bool compareToCPU = false;

bool QATest = false;

int blockSize = 256;

bool useHostMem = false;

bool fp64 = false;

bool useCpu = false;

int numDevsRequested = 1;

bool displayEnabled = true;

bool bPause = false;

bool bFullscreen = false;

bool bDispInteractions = false;

bool bSupportDouble = false;

int flopsPerInteraction = 20;

char deviceName[100];

enum { M_VIEW = 0, M_MOVE };

int numBodies = 16384;

std::string tipsyFile = "";

int numIterations = 0; // run until exit

void computePerfStats(double &interactionsPerSecond, double &gflops,

float milliseconds, int iterations)

{

// Dubla precizie foloseste funcționarea intrinsecă urmată de rafinament,
     // Nr.de functionare;
     // (Notă Astrofizicienii folosesc 38 flops pe interacțiune)
     // Bazat pe "precedent istoric", dar ei folosesc FLOP / s ;
     // Măsură de "transfer știință". Suntem folosind-o ca o măsură de
     // Transfer hardware.
     // Const int flopsPerInteraction = fp64? 30: 20; interactionsPerSecond = (float)numBodies * (float)numBodies;

interactionsPerSecond *= 1e-9 * iterations * 1000 / milliseconds;

gflops = interactionsPerSecond * (float)flopsPerInteraction;

}

////////////////////////////////////////

// Demo Parametrii

////////////////////////////////////////

struct NBodyParams

{

float m_timestep;

float m_clusterScale;

float m_velocityScale;

float m_softening;

float m_damping;

float m_pointSize;

float m_x, m_y, m_z;

void print()

{

printf("{ %f, %f, %f, %f, %f, %f, %f, %f, %f },\n",

m_timestep, m_clusterScale, m_velocityScale,

m_softening, m_damping, m_pointSize, m_x, m_y, m_z);

}

};

NBodyParams demoParams[] =

{

{ 0.016f, 1.54f, 8.0f, 0.1f, 1.0f, 1.0f, 0, -2, -100},

{ 0.016f, 0.68f, 20.0f, 0.1f, 1.0f, 0.8f, 0, -2, -30},

{ 0.0006f, 0.16f, 1000.0f, 1.0f, 1.0f, 0.07f, 0, 0, -1.5f},

{ 0.0006f, 0.16f, 1000.0f, 1.0f, 1.0f, 0.07f, 0, 0, -1.5f},

{ 0.0019f, 0.32f, 276.0f, 1.0f, 1.0f, 0.07f, 0, 0, -5},

{ 0.0016f, 0.32f, 272.0f, 0.145f, 1.0f, 0.08f, 0, 0, -5},

{ 0.016000f, 6.040000f, 0.000000f, 1.000000f, 1.000000f, 0.760000f, 0, 0, -50},

};

int numDemos = sizeof(demoParams) / sizeof(NBodyParams);

bool cycleDemo = true;

int activeDemo = 0;

float demoTime = 10000.0f; // ms

StopWatchInterface *demoTimer = NULL, *timer = NULL;

// executa mai multe iteratii pentru a folosi un timp mediu

NBodyParams activeParams = demoParams[activeDemo];

// The UI.

ParamListGL *paramlist; // lista de parametrii

bool bShowSliders = true;

// fps

static int fpsCount = 0;

static int fpsLimit = 5;

cudaEvent_t startEvent, stopEvent;

cudaEvent_t hostMemSyncEvent;

template <typename T>

class NBodyDemo

{

public:

static void Create()

{

m_singleton = new NBodyDemo;

}

static void Destroy()

{

delete m_singleton;

}

static void init(int numBodies, int numDevices, int blockSize,

bool usePBO, bool useHostMem, bool useCpu)

{

m_singleton->_init(numBodies, numDevices, blockSize, usePBO, useHostMem, useCpu);

}

static void reset(int numBodies, NBodyConfig config)

{

m_singleton->_reset(numBodies, config);

}

static void selectDemo(int index)

{

m_singleton->_selectDemo(index);

}

static bool compareResults(int numBodies)

{

return m_singleton->_compareResults(numBodies);

}

static void runBenchmark(int iterations)

{

m_singleton->_runBenchmark(iterations);

}

static void updateParams()

{

m_singleton->m_nbody->setSoftening(activeParams.m_softening);

m_singleton->m_nbody->setDamping(activeParams.m_damping);

}

static void updateSimulation()

{

m_singleton->m_nbody->update(activeParams.m_timestep);

}

static void display()

{

m_singleton->m_renderer->setSpriteSize(activeParams.m_pointSize);

if (useHostMem)

{

// Aceasta sincronizare eveniment este necesară pentru că redam din memoria gazda
                 // Scris. Dacă nu vom aștepta până CUDA se face actualizarea, vom face o actualizare parțiala // update de date

if (!useCpu)

{

cudaEventSynchronize(hostMemSyncEvent);

}

m_singleton->m_renderer->setPositions(

m_singleton->m_nbody->getArray(BODYSYSTEM_POSITION),

m_singleton->m_nbody->getNumBodies());

}

else

{

m_singleton->m_renderer->setPBO(m_singleton->m_nbody->getCurrentReadBuffer(),

m_singleton->m_nbody->getNumBodies(),

(sizeof(T) > 4));

}

// afisare particule

m_singleton->m_renderer->display(displayMode);

}

static void getArrays(T *pos, T *vel)

{

T *_pos = m_singleton->m_nbody->getArray(BODYSYSTEM_POSITION);

T *_vel = m_singleton->m_nbody->getArray(BODYSYSTEM_VELOCITY);

memcpy(pos, _pos, m_singleton->m_nbody->getNumBodies() * 4 * sizeof(T));

memcpy(vel, _vel, m_singleton->m_nbody->getNumBodies() * 4 * sizeof(T));

}

static void setArrays(const T *pos, const T *vel)

{

if (pos != m_singleton->m_hPos)

{

memcpy(m_singleton->m_hPos, pos, numBodies * 4 * sizeof(T));

}

if (vel != m_singleton->m_hVel)

{

memcpy(m_singleton->m_hVel, vel, numBodies * 4 * sizeof(T));

}

m_singleton->m_nbody->setArray(BODYSYSTEM_POSITION, m_singleton->m_hPos);

m_singleton->m_nbody->setArray(BODYSYSTEM_VELOCITY, m_singleton->m_hVel);

if (!benchmark && !useCpu && !compareToCPU)

{

m_singleton->_resetRenderer();

}

}

private:

static NBodyDemo *m_singleton;

BodySystem<T> *m_nbody;

BodySystemCUDA<T> *m_nbodyCuda;

BodySystemCPU<T> *m_nbodyCpu;

ParticleRenderer *m_renderer;

T *m_hPos;

T *m_hVel;

float *m_hColor;

private:

NBodyDemo()

: m_nbody(0),

m_nbodyCuda(0),

m_nbodyCpu(0),

m_renderer(0),

m_hPos(0),

m_hVel(0),

m_hColor(0)

{

}

~NBodyDemo()

{

if (m_nbodyCpu)

{

delete m_nbodyCpu;

}

if (m_nbodyCuda)

{

delete m_nbodyCuda;

}

if (m_hPos)

{

delete [] m_hPos;

}

if (m_hVel)

{

delete [] m_hVel;

}

if (m_hColor)

{

delete [] m_hColor;

}

sdkDeleteTimer(&demoTimer);

if (!benchmark && !compareToCPU)

delete m_renderer;

}

void _init(int numBodies, int numDevices, int blockSize,

bool bUsePBO, bool useHostMem, bool useCpu)

{

if (useCpu)

{

m_nbodyCpu = new BodySystemCPU<T>(numBodies);

m_nbody = m_nbodyCpu;

m_nbodyCuda = 0;

}

else

{

m_nbodyCuda = new BodySystemCUDA<T>(numBodies, numDevices, blockSize, bUsePBO, useHostMem);

m_nbody = m_nbodyCuda;

m_nbodyCpu = 0;

}

// alocare memorie gazda

m_hPos = new T[numBodies*4];

m_hVel = new T[numBodies*4];

m_hColor = new float[numBodies*4];

m_nbody->setSoftening(activeParams.m_softening);

m_nbody->setDamping(activeParams.m_damping);

if (useCpu)

{

sdkCreateTimer(&timer);

sdkStartTimer(&timer);

}

else

{

checkCudaErrors(cudaEventCreate(&startEvent));

checkCudaErrors(cudaEventCreate(&stopEvent));

checkCudaErrors(cudaEventCreate(&hostMemSyncEvent));

}

if (!benchmark && !compareToCPU)

{

m_renderer = new ParticleRenderer;

_resetRenderer();

}

sdkCreateTimer(&demoTimer);

sdkStartTimer(&demoTimer);

}

void _reset(int numBodies, NBodyConfig config)

{

if (tipsyFile == "")

{

randomizeBodies(config, m_hPos, m_hVel, m_hColor,

activeParams.m_clusterScale,

activeParams.m_velocityScale,

numBodies, true);

setArrays(m_hPos, m_hVel);

}

else

{

m_nbody->loadTipsyFile(tipsyFile);

::numBodies = m_nbody->getNumBodies();

}

}

void _resetRenderer()

{

if (fp64)

{

float color[4] = { 0.4f, 0.8f, 0.1f, 1.0f};

m_renderer->setBaseColor(color);

}

else

{

float color[4] = { 1.0f, 0.6f, 0.3f, 1.0f};

m_renderer->setBaseColor(color);

}

m_renderer->setColors(m_hColor, m_nbody->getNumBodies());

m_renderer->setSpriteSize(activeParams.m_pointSize);

}

void _selectDemo(int index)

{

assert(index < numDemos);

activeParams = demoParams[index];

camera_trans[0] = camera_trans_lag[0] = activeParams.m_x;

camera_trans[1] = camera_trans_lag[1] = activeParams.m_y;

camera_trans[2] = camera_trans_lag[2] = activeParams.m_z;

reset(numBodies, NBODY_CONFIG_SHELL);

sdkResetTimer(&demoTimer);

}

bool _compareResults(int numBodies)

{

assert(m_nbodyCuda);

bool passed = true;

m_nbody->update(0.001f);

{

m_nbodyCpu = new BodySystemCPU<T>(numBodies);

m_nbodyCpu->setArray(BODYSYSTEM_POSITION, m_hPos);

m_nbodyCpu->setArray(BODYSYSTEM_VELOCITY, m_hVel);

m_nbodyCpu->update(0.001f);

T *cudaPos = m_nbodyCuda->getArray(BODYSYSTEM_POSITION);

T *cpuPos = m_nbodyCpu->getArray(BODYSYSTEM_POSITION);

T tolerance = 0.0005f;

for (int i = 0; i < numBodies; i++)

{

if (fabs(cpuPos[i] – cudaPos[i]) > tolerance)

{

passed = false;

printf("Error: (host)%f != (device)%f\n", cpuPos[i], cudaPos[i]);

}

}

}

return passed;

}

void _runBenchmark(int iterations)

{

if (!useCpu)

{

m_nbody->update(activeParams.m_timestep);

}

if (useCpu)

{

sdkCreateTimer(&timer);

sdkStartTimer(&timer);

}

else

{

checkCudaErrors(cudaEventRecord(startEvent, 0));

}

for (int i = 0; i < iterations; ++i)

{

m_nbody->update(activeParams.m_timestep);

}

float milliseconds = 0;

if (useCpu)

{

sdkStopTimer(&timer);

milliseconds = sdkGetTimerValue(&timer);

sdkStartTimer(&timer);

}

else

{

checkCudaErrors(cudaEventRecord(stopEvent, 0));

checkCudaErrors(cudaEventSynchronize(stopEvent));

checkCudaErrors(cudaEventElapsedTime(&milliseconds, startEvent, stopEvent));

}

double interactionsPerSecond = 0;

double gflops = 0;

computePerfStats(interactionsPerSecond, gflops, milliseconds, iterations);

printf("%d bodies, total time for %d iterations: %.3f ms\n",

numBodies, iterations, milliseconds);

printf("= %.3f billion interactions per second\n", interactionsPerSecond);

printf("= %.3f %s-precision GFLOP/s at %d flops per interaction\n", gflops,

(sizeof(T) > 4) ? "double" : "single", flopsPerInteraction);

}

};

void finalize()

{

if (!useCpu)

{

checkCudaErrors(cudaEventDestroy(startEvent));

checkCudaErrors(cudaEventDestroy(stopEvent));

checkCudaErrors(cudaEventDestroy(hostMemSyncEvent));

}

NBodyDemo<float>::Destroy();

if (bSupportDouble) NBodyDemo<double>::Destroy();

}

template <> NBodyDemo<double> *NBodyDemo<double>::m_singleton = 0;

template <> NBodyDemo<float> *NBodyDemo<float>::m_singleton = 0;

template <typename T_new, typename T_old>

void switchDemoPrecision()

{

cudaDeviceSynchronize();

fp64 = !fp64;

flopsPerInteraction = fp64 ? 30 : 20;

T_old *oldPos = new T_old[numBodies * 4];

T_old *oldVel = new T_old[numBodies * 4];

NBodyDemo<T_old>::getArrays(oldPos, oldVel);

// convertire float in double

T_new *newPos = new T_new[numBodies * 4];

T_new *newVel = new T_new[numBodies * 4];

for (int i = 0; i < numBodies * 4; i++)

{

newPos[i] = (T_new)oldPos[i];

newVel[i] = (T_new)oldVel[i];

}

NBodyDemo<T_new>::setArrays(newPos, newVel);

cudaDeviceSynchronize();

delete [] oldPos;

delete [] oldVel;

delete [] newPos;

delete [] newVel;

}

// verifica erori pt OpenGL

inline

void checkGLErrors(const char *s)

{

GLenum error;

while ((error = glGetError()) != GL_NO_ERROR)

{

fprintf(stderr, "%s: error – %s\n", s, (char *) gluErrorString(error));

}

}

void initGL(int *argc, char **argv)

{

// performanta CUDA – folosind OPENGL

glutInit(argc, argv);

glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE);

glutInitWindowSize(720, 480);

glutCreateWindow("CUDA n-body system");

if (bFullscreen)

{

glutFullScreen();

}

GLenum err = glewInit();

if (GLEW_OK != err)

{

printf("GLEW Error: %s\n", glewGetErrorString(err));

// curatare inainte de a iesi din aplicatie

cudaDeviceReset();

exit(EXIT_FAILURE);

}

else if (!glewIsSupported("GL_VERSION_2_0 "

"GL_VERSION_1_5 "

"GL_ARB_multitexture "

"GL_ARB_vertex_buffer_object"))

{

fprintf(stderr, "Required OpenGL extensions missing.");

exit(EXIT_FAILURE);

}

else

{

#if defined(WIN32)

wglSwapIntervalEXT(0);

#elif defined(LINUX)

glxSwapIntervalSGI(0);

#endif

}

glEnable(GL_DEPTH_TEST);

glClearColor(0.0, 0.0, 0.0, 1.0);

checkGLErrors("initGL");

}

void initParameters()

{

// creare o noua lista de parametrii

paramlist = new ParamListGL("sliders");

paramlist->SetBarColorInner(0.8f, 0.8f, 0.0f);

// Adaugari anumiti parametrii la lista

// Dimensiune punct

paramlist->AddParam(new Param<float>("Point Size", activeParams.m_pointSize,

0.001f, 10.0f, 0.01f, &activeParams.m_pointSize));

// Viteza de amortizare

paramlist->AddParam(new Param<float>("Velocity Damping", activeParams.m_damping,

0.5f, 1.0f, .0001f, &(activeParams.m_damping)));

// Factorul de dedurizare

paramlist->AddParam(new Param<float>("Softening Factor", activeParams.m_softening,

0.001f, 1.0f, .0001f, &(activeParams.m_softening)));

// Mărimea pasului de timp paramlist->AddParam(new Param<float>("Time Step", activeParams.m_timestep,

0.0f, 1.0f, .0001f, &(activeParams.m_timestep)));

// SCALA CLUSTER

paramlist->AddParam(new Param<float>("Cluster Scale", activeParams.m_clusterScale,

0.0f, 10.0f, 0.01f, &(activeParams.m_clusterScale)));

// Scala vitezei

paramlist->AddParam(new Param<float>("Velocity Scale", activeParams.m_velocityScale,

0.0f, 1000.0f, 0.1f, &activeParams.m_velocityScale));

}

void selectDemo(int activeDemo)

{

if (fp64)

{

NBodyDemo<double>::selectDemo(activeDemo);

}

else

{

NBodyDemo<float>::selectDemo(activeDemo);

}

}

void updateSimulation()

{

if (fp64)

{

NBodyDemo<double>::updateSimulation();

}

else

{

NBodyDemo<float>::updateSimulation();

}

}

void displayNBodySystem()

{

if (fp64)

{

NBodyDemo<double>::display();

}

else

{

NBodyDemo<float>::display();

}

}

void display()

{

static double gflops = 0;

static double ifps = 0;

static double interactionsPerSecond = 0;

// updatare

if (!bPause)

{

if (cycleDemo && (sdkGetTimerValue(&demoTimer) > demoTime))

{

activeDemo = (activeDemo + 1) % numDemos;

selectDemo(activeDemo);

}

updateSimulation();

if (!useCpu)

{

cudaEventRecord(hostMemSyncEvent, 0);

}

}

glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

if (displayEnabled)

{

// Vizualizare transformare

{

glMatrixMode(GL_MODELVIEW);

glLoadIdentity();

for (int c = 0; c < 3; ++c)

{

camera_trans_lag[c] += (camera_trans[c] – camera_trans_lag[c]) * inertia;

camera_rot_lag[c] += (camera_rot[c] – camera_rot_lag[c]) * inertia;

}

glTranslatef(camera_trans_lag[0], camera_trans_lag[1], camera_trans_lag[2]);

glRotatef(camera_rot_lag[0], 1.0, 0.0, 0.0);

glRotatef(camera_rot_lag[1], 0.0, 1.0, 0.0);

}

displayNBodySystem();

// afișare interfață cu utilizatorul

if (bShowSliders)

{

glBlendFunc(GL_ONE_MINUS_DST_COLOR, GL_ZERO); // inversare culori

glEnable(GL_BLEND);

paramlist->Render(0, 0);

glDisable(GL_BLEND);

}

if (bFullscreen)

{

beginWinCoords();

char msg0[256], msg1[256], msg2[256];

if (bDispInteractions)

{

sprintf(msg1, "%0.2f billion interactions per second", interactionsPerSecond);

}

else

{

sprintf(msg1, "%0.2f GFLOP/s", gflops);

}

sprintf(msg0, "%s", deviceName);

sprintf(msg2, "%0.2f FPS [%s | %d bodies]",

ifps, fp64 ? "double precision" : "single precision", numBodies);

glBlendFunc(GL_ONE_MINUS_DST_COLOR, GL_ZERO); // inversare culori

glEnable(GL_BLEND);

glColor3f(0.46f, 0.73f, 0.0f);

glPrint(80, glutGet(GLUT_WINDOW_HEIGHT) – 122, msg0, GLUT_BITMAP_TIMES_ROMAN_24);

glColor3f(1.0f, 1.0f, 1.0f);

glPrint(80, glutGet(GLUT_WINDOW_HEIGHT) – 96, msg2, GLUT_BITMAP_TIMES_ROMAN_24);

glColor3f(1.0f, 1.0f, 1.0f);

glPrint(80, glutGet(GLUT_WINDOW_HEIGHT) – 70, msg1, GLUT_BITMAP_TIMES_ROMAN_24);

glDisable(GL_BLEND);

endWinCoords();

}

glutSwapBuffers();

}

fpsCount++;

// rata de frame

if (fpsCount >= fpsLimit)

{

char fps[256];

float milliseconds = 1;

// oprire timp

if (useCpu)

{

milliseconds = sdkGetTimerValue(&timer);

sdkResetTimer(&timer);

}

else

{

checkCudaErrors(cudaEventRecord(stopEvent, 0));

checkCudaErrors(cudaEventSynchronize(stopEvent));

checkCudaErrors(cudaEventElapsedTime(&milliseconds, startEvent, stopEvent));

}

milliseconds /= (float)fpsCount;

computePerfStats(interactionsPerSecond, gflops, milliseconds, 1);

ifps = 1.f / (milliseconds / 1000.f);

sprintf(fps,

"CUDA N-Body (%d bodies): "

"%0.1f fps | %0.1f BIPS | %0.1f GFLOP/s | %s",

numBodies, ifps, interactionsPerSecond, gflops,

fp64 ? "double precision" : "single precision");

glutSetWindowTitle(fps);

fpsCount = 0;

fpsLimit = (ifps > 1.f) ? (int)ifps : 1;

if (bPause)

{

fpsLimit = 0;

}

// restart timp

if (!useCpu)

{

checkCudaErrors(cudaEventRecord(startEvent, 0));

}

}

glutReportErrors();

}

void reshape(int w, int h)

{

glMatrixMode(GL_PROJECTION);

glLoadIdentity();

gluPerspective(60.0, (float) w / (float) h, 0.1, 1000.0);

glMatrixMode(GL_MODELVIEW);

glViewport(0, 0, w, h);

}

void updateParams()

{

if (fp64)

{

NBodyDemo<double>::updateParams();

}

else

{

NBodyDemo<float>::updateParams();

}

}

void mouse(int button, int state, int x, int y)

{

if (bShowSliders)

{

// functia de apel – mouse

if (paramlist->Mouse(x, y, button, state))

{

updateParams();

}

}

int mods;

if (state == GLUT_DOWN)

{

buttonState |= 1<<button;

}

else if (state == GLUT_UP)

{

buttonState = 0;

}

mods = glutGetModifiers();

if (mods & GLUT_ACTIVE_SHIFT)

{

buttonState = 2;

}

else if (mods & GLUT_ACTIVE_CTRL)

{

buttonState = 3;

}

ox = x;

oy = y;

glutPostRedisplay();

}

void motion(int x, int y)

{

if (bShowSliders)

{

// apel la parametrii

if (paramlist->Motion(x, y))

{

updateParams();

glutPostRedisplay();

return;

}

}

float dx = (float)(x – ox);

float dy = (float)(y – oy);

if (buttonState == 3)

{

// stanga + mijloc = zoom

camera_trans[2] += (dy / 100.0f) * 0.5f * fabs(camera_trans[2]);

}

else if (buttonState & 2)

{

// mijloc – traduce

camera_trans[0] += dx / 100.0f;

camera_trans[1] -= dy / 100.0f;

}

else if (buttonState & 1)

{

// stanga – roteste

camera_rot[0] += dy / 5.0f;

camera_rot[1] += dx / 5.0f;

}

ox = x;

oy = y;

glutPostRedisplay();

}

void key(unsigned char key, int /*x*/, int /*y*/)

{

switch (key)

{

case ' ':

bPause = !bPause;

break;

case 27: // escape

case 'q':

case 'Q':

finalize();

//finalizare -> reset

cudaDeviceReset();

exit(EXIT_SUCCESS);

break;

case 13: // return

if (bSupportDouble)

{

if (fp64)

{

switchDemoPrecision<float, double>();

}

else

{

switchDemoPrecision<double, float>();

}

printf("> %s precision floating point simulation\n", fp64 ? "Double" : "Single");

}

break;

case '`':

bShowSliders = !bShowSliders;

break;

case 'g':

case 'G':

bDispInteractions = !bDispInteractions;

break;

case 'p':

case 'P':

displayMode = (ParticleRenderer::DisplayMode)

((displayMode + 1) % ParticleRenderer::PARTICLE_NUM_MODES);

break;

case 'c':

case 'C':

cycleDemo = !cycleDemo;

printf("Cycle Demo Parameters: %s\n", cycleDemo ? "ON" : "OFF");

break;

case '[':

activeDemo = (activeDemo == 0) ? numDemos – 1 : (activeDemo – 1) % numDemos;

selectDemo(activeDemo);

break;

case ']':

activeDemo = (activeDemo + 1) % numDemos;

selectDemo(activeDemo);

break;

case 'd':

case 'D':

displayEnabled = !displayEnabled;

break;

case 'o':

case 'O':

activeParams.print();

break;

case '1':

if (fp64)

{

NBodyDemo<double>::reset(numBodies, NBODY_CONFIG_SHELL);

}

else

{

NBodyDemo<float>::reset(numBodies, NBODY_CONFIG_SHELL);

}

break;

case '2':

if (fp64)

{

NBodyDemo<double>::reset(numBodies, NBODY_CONFIG_RANDOM);

}

else

{

NBodyDemo<float>::reset(numBodies, NBODY_CONFIG_RANDOM);

}

break;

case '3':

if (fp64)

{

NBodyDemo<double>::reset(numBodies, NBODY_CONFIG_EXPAND);

}

else

{

NBodyDemo<float>::reset(numBodies, NBODY_CONFIG_EXPAND);

}

break;

}

glutPostRedisplay();

}

void special(int key, int x, int y)

{

paramlist->Special(key, x, y);

glutPostRedisplay();

}

void idle(void)

{

glutPostRedisplay();

}

void

showHelp()

{

printf("\t-fullscreen (run n-body simulation in fullscreen mode)\n");

printf("\t-fp64 (use double precision floating point values for simulation)\n");

printf("\t-hostmem (stores simulation data in host memory)\n");

printf("\t-benchmark (run benchmark to measure performance) \n");

printf("\t-numbodies=<N> (number of bodies (>= 1) to run in simulation) \n");

printf("\t-device=<d> (where d=0,1,2…. for the CUDA device to use)\n");

printf("\t-numdevices=<i> (where i=(number of CUDA devices > 0) to use for simulation)\n");

printf("\t-compare (compares simulation results running once on the default GPU and once on the CPU)\n");

printf("\t-cpu (run n-body simulation on the CPU)\n");

printf("\t-tipsy=<file.bin> (load a tipsy model file for simulation)\n\n");

}

//////////////////////////////////////////////////////////////////////////////

// MAIN PROGRAM

//////////////////////////////////////////////////////////////////////////////

int

main(int argc, char **argv)

{

bool bTestResults = true;

#if defined(__linux__)

setenv ("DISPLAY", ":0", 0);

#endif

if (checkCmdLineFlag(argc, (const char **)argv, "help"))

{

printf("\n> Command line options\n");

showHelp();

return 0;

}

printf("Run \"nbody -benchmark [-numbodies=<numBodies>]\" to measure perfomance.\n");

showHelp();

printf("NOTE: The CUDA Samples are not meant for performance measurements. Results may vary when GPU Boost is enabled.\n\n");

bFullscreen = (checkCmdLineFlag(argc, (const char **) argv, "fullscreen") != 0);

if (bFullscreen)

{

bShowSliders = false;

}

benchmark = (checkCmdLineFlag(argc, (const char **) argv, "benchmark") != 0);

compareToCPU = ((checkCmdLineFlag(argc, (const char **) argv, "compare") != 0) ||

(checkCmdLineFlag(argc, (const char **) argv, "qatest") != 0));

QATest = (checkCmdLineFlag(argc, (const char **) argv, "qatest") != 0);

useHostMem = (checkCmdLineFlag(argc, (const char **) argv, "hostmem") != 0);

fp64 = (checkCmdLineFlag(argc, (const char **) argv, "fp64") != 0);

flopsPerInteraction = fp64 ? 30 : 20;

useCpu = (checkCmdLineFlag(argc, (const char **) argv, "cpu") != 0);

if (checkCmdLineFlag(argc, (const char **)argv, "numdevices"))

{

numDevsRequested = getCmdLineArgumentInt(argc, (const char **) argv, "numdevices");

if (numDevsRequested < 1)

{

printf("Error: \"number of CUDA devices\" specified %d is invalid. Value should be >= 1\n", numDevsRequested);

exit(bTestResults ? EXIT_SUCCESS : EXIT_FAILURE);

}

else

{

printf("number of CUDA devices = %d\n", numDevsRequested);

}

}

// date prin memoria gazda

if (numDevsRequested > 1)

{

useHostMem = true;

}

int numDevsAvailable = 0;

bool customGPU = false;

cudaGetDeviceCount(&numDevsAvailable);

if (numDevsAvailable < numDevsRequested)

{

printf("Error: only %d Devices available, %d requested. Exiting.\n", numDevsAvailable, numDevsRequested);

exit(EXIT_SUCCESS);

}

printf("> %s mode\n", bFullscreen ? "Fullscreen" : "Windowed");

printf("> Simulation data stored in %s memory\n", useHostMem ? "system" : "video");

printf("> %s precision floating point simulation\n", fp64 ? "Double" : "Single");

printf("> %d Devices used for simulation\n", numDevsRequested);

int devID;

cudaDeviceProp props;

if (useCpu)

{

useHostMem = true;

compareToCPU = false;

bSupportDouble = true;

#ifdef OPENMP

printf("> Simulation with CPU using OpenMP\n");

#else

printf("> Simulation with CPU\n");

#endif

}

// INITIALIZARE LIBRARIA GLUT daca este necesar

if (!benchmark && !compareToCPU)

{

initGL(&argc, argv);

initParameters();

}

if (!useCpu)

{

if (benchmark || compareToCPU || useHostMem)

{

if (checkCmdLineFlag(argc, (const char **)argv, "device"))

{

customGPU = true;

}

devID = findCudaDevice(argc, (const char **)argv);

}

else // or with GL interop:

{

if (checkCmdLineFlag(argc, (const char **)argv, "device"))

{

customGPU = true;

}

devID = findCudaGLDevice(argc, (const char **)argv);

}

checkCudaErrors(cudaGetDevice(&devID));

checkCudaErrors(cudaGetDeviceProperties(&props, devID));

bSupportDouble = true;

#if CUDART_VERSION < 4000

if (numDevsRequested > 1)

{

printf("MultiGPU n-body requires CUDA 4.0 or later\n");

// curatare

cudaDeviceReset();

exit(EXIT_SUCCESS);

}

#endif

// Initializare DEVICE

if (numDevsRequested > 1 && customGPU)

{

printf("You can't use –numdevices and –device at the same time.\n");

exit(EXIT_SUCCESS);

}

if (customGPU || numDevsRequested == 1)

{

cudaDeviceProp props;

checkCudaErrors(cudaGetDeviceProperties(&props, devID));

printf("> Compute %d.%d CUDA device: [%s]\n", props.major, props.minor, props.name);

}

else

{

for (int i = 0; i < numDevsRequested; i++)

{

cudaDeviceProp props;

checkCudaErrors(cudaGetDeviceProperties(&props, i));

printf("> Compute %d.%d CUDA device: [%s]\n", props.major, props.minor, props.name);

if (useHostMem)

{

#if CUDART_VERSION >= 2020

if (!props.canMapHostMemory)

{

fprintf(stderr, "Device %d cannot map host memory!\n", devID);

cudaDeviceReset();

exit(EXIT_SUCCESS);

}

if (numDevsRequested > 1)

{

checkCudaErrors(cudaSetDevice(i));

}

checkCudaErrors(cudaSetDeviceFlags(cudaDeviceMapHost));

#else

fprintf(stderr, "This CUDART version does not support <cudaDeviceProp.canMapHostMemory> field\n");

cudaDeviceReset();

exit(EXIT_SUCCESS);

#endif

}

}

if (props.major*10 + props.minor <= 12)

{

bSupportDouble = false;

}

}

//if(numDevsRequested > 1)

// checkCudaErrors(cudaSetDevice(devID));

if (fp64 && !bSupportDouble)

{

fprintf(stderr, "One or more of the requested devices does not support double precision floating-point\n");

cudaDeviceReset();

exit(EXIT_SUCCESS);

}

}

numIterations = 0;

blockSize = 0;

if (checkCmdLineFlag(argc, (const char **)argv, "i"))

{

numIterations = getCmdLineArgumentInt(argc, (const char **)argv, "i");

}

if (checkCmdLineFlag(argc, (const char **) argv, "blockSize"))

{

blockSize = getCmdLineArgumentInt(argc, (const char **)argv, "blockSize");

}

if (blockSize == 0) // dimensiunea blocului

blockSize = 256;

NR CORPURI

if (useCpu)

#ifdef OPENMP

numBodies = 8192;

#else

numBodies = 4096;

#endif

else if (numDevsRequested == 1)

{

numBodies = compareToCPU ? 4096 : blockSize*4*props.multiProcessorCount;

}

else

{

numBodies = 0;

for (int i = 0; i < numDevsRequested; i++)

{

cudaDeviceProp props;

checkCudaErrors(cudaGetDeviceProperties(&props, i));

numBodies += blockSize*(props.major >= 2 ? 4 : 1)*props.multiProcessorCount;

}

}

if (checkCmdLineFlag(argc, (const char **) argv, "numbodies"))

{

numBodies = getCmdLineArgumentInt(argc, (const char **)argv, "numbodies");

if (numBodies < 1)

{

printf("Error: \"number of bodies\" specified %d is invalid. Value should be >= 1\n", numBodies);

exit(bTestResults ? EXIT_SUCCESS : EXIT_FAILURE);

}

else if (numBodies % blockSize)

{

int newNumBodies = ((numBodies / blockSize) + 1) * blockSize;

printf("Warning: \"number of bodies\" specified %d is not a multiple of %d.\n", numBodies, blockSize);

printf("Rounding up to the nearest multiple: %d.\n", newNumBodies);

numBodies = newNumBodies;

}

else

{

printf("number of bodies = %d\n", numBodies);

}

}

char *fname;

if (getCmdLineArgumentString(argc, (const char **)argv, "tipsy", &fname))

{

tipsyFile.assign(fname, strlen(fname));

cycleDemo = false;

bShowSliders = false;

}

if (numBodies <= 1024)

{

activeParams.m_clusterScale = 1.52f;

activeParams.m_velocityScale = 2.f;

}

else if (numBodies <= 2048)

{

activeParams.m_clusterScale = 1.56f;

activeParams.m_velocityScale = 2.64f;

}

else if (numBodies <= 4096)

{

activeParams.m_clusterScale = 1.68f;

activeParams.m_velocityScale = 2.98f;

}

else if (numBodies <= 8192)

{

activeParams.m_clusterScale = 1.98f;

activeParams.m_velocityScale = 2.9f;

}

else if (numBodies <= 16384)

{

activeParams.m_clusterScale = 1.54f;

activeParams.m_velocityScale = 8.f;

}

else if (numBodies <= 32768)

{

activeParams.m_clusterScale = 1.44f;

activeParams.m_velocityScale = 11.f;

}

// Creare demo si punere in aplicare

NBodyDemo<float>::Create();

NBodyDemo<float>::init(numBodies, numDevsRequested, blockSize, !(benchmark || compareToCPU || useHostMem), useHostMem, useCpu);

NBodyDemo<float>::reset(numBodies, NBODY_CONFIG_SHELL);

if (bSupportDouble)

{

NBodyDemo<double>::Create();

NBodyDemo<double>::init(numBodies, numDevsRequested, blockSize, !(benchmark || compareToCPU || useHostMem), useHostMem, useCpu);

NBodyDemo<double>::reset(numBodies, NBODY_CONFIG_SHELL);

}

if (fp64)

{

if (benchmark)

{

if (numIterations <= 0)

{

numIterations = 10;

}

else if (numIterations > 10)

{

printf("Advisory: setting a high number of iterations\n");

printf("in benchmark mode may cause failure on Windows\n");

printf("Vista and Win7. On these OSes, set iterations <= 10\n");

}

NBodyDemo<double>::runBenchmark(numIterations);

}

else if (compareToCPU)

{

bTestResults = NBodyDemo<double>::compareResults(numBodies);

}

else

{

glutDisplayFunc(display);

glutReshapeFunc(reshape);

glutMouseFunc(mouse);

glutMotionFunc(motion);

glutKeyboardFunc(key);

glutSpecialFunc(special);

glutIdleFunc(idle);

if (!useCpu)

{

checkCudaErrors(cudaEventRecord(startEvent, 0));

}

glutMainLoop();

}

}

else

{

if (benchmark)

{

if (numIterations <= 0)

{

numIterations = 10;

}

NBodyDemo<float>::runBenchmark(numIterations);

}

else if (compareToCPU)

{

bTestResults = NBodyDemo<float>::compareResults(numBodies);

}

else

{

glutDisplayFunc(display);

glutReshapeFunc(reshape);

glutMouseFunc(mouse);

glutMotionFunc(motion);

glutKeyboardFunc(key);

glutSpecialFunc(special);

glutIdleFunc(idle);

if (!useCpu)

{

checkCudaErrors(cudaEventRecord(startEvent, 0));

}

glutMainLoop();

}

}

finalize();

// curatare inainte de iesirea din aplicatie

cudaDeviceReset();

exit(bTestResults ? EXIT_SUCCESS : EXIT_FAILURE);

}

“render_particles.cpp”

#include "render_particles.h"

#include <GL/glew.h>

#include <cuda_runtime.h>

#include <cuda_gl_interop.h>

#include <helper_cuda.h>

#include <helper_cuda_gl.h>

#if defined(__APPLE__) || defined(MACOSX)

#pragma clang diagnostic ignored "-Wdeprecated-declarations"

#include <GLUT/glut.h>

#else

#include <GL/freeglut.h>

#endif

#include <math.h>

#include <assert.h>

#define GL_POINT_SPRITE_ARB 0x8861

#define GL_COORD_REPLACE_ARB 0x8862

#define GL_VERTEX_PROGRAM_POINT_SIZE_NV 0x8642

ParticleRenderer::ParticleRenderer()

: m_pos(0),

m_numParticles(0),

m_pointSize(1.0f),

m_spriteSize(2.0f),

m_vertexShader(0),

m_vertexShaderPoints(0),

m_pixelShader(0),

m_programPoints(0),

m_programSprites(0),

m_texture(0),

m_pbo(0),

m_vboColor(0),

m_bFp64Positions(false)

{

_initGL();

}

ParticleRenderer::~ParticleRenderer()

{

m_pos = 0;

}

void ParticleRenderer::resetPBO()

{

glDeleteBuffers(1, (GLuint *)&m_pbo);

}

void ParticleRenderer::setPositions(float *pos, int numParticles)

{

m_pos = pos;

m_numParticles = numParticles;

if (!m_pbo)

{

glGenBuffers(1, (GLuint *)&m_pbo);

}

glBindBuffer(GL_ARRAY_BUFFER_ARB, m_pbo);

glBufferData(GL_ARRAY_BUFFER_ARB, numParticles * 4 * sizeof(float), pos, GL_STATIC_DRAW_ARB);

glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

SDK_CHECK_ERROR_GL();

}

void ParticleRenderer::setPositions(double *pos, int numParticles)

{

m_bFp64Positions = true;

m_pos_fp64 = pos;

m_numParticles = numParticles;

if (!m_pbo)

{

glGenBuffers(1, (GLuint *)&m_pbo);

}

glBindBuffer(GL_ARRAY_BUFFER_ARB, m_pbo);

glBufferData(GL_ARRAY_BUFFER_ARB, numParticles * 4 * sizeof(double), pos, GL_STATIC_DRAW_ARB);

glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

SDK_CHECK_ERROR_GL();

}

void ParticleRenderer::setColors(float *color, int numParticles)

{

glBindBuffer(GL_ARRAY_BUFFER_ARB, m_vboColor);

glBufferData(GL_ARRAY_BUFFER_ARB, numParticles * 4 * sizeof(float), color, GL_STATIC_DRAW_ARB);

glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

}

void ParticleRenderer::setBaseColor(float color[4])

{

for (int i = 0; i < 4; i++)

m_baseColor[i] = color[i];

}

void ParticleRenderer::setPBO(unsigned int pbo, int numParticles, bool fp64)

{

m_pbo = pbo;

m_numParticles = numParticles;

if (fp64)

m_bFp64Positions = true;

}

void ParticleRenderer::_drawPoints(bool color)

{

if (!m_pbo)

{

glBegin(GL_POINTS);

{

int k = 0;

for (int i = 0; i < m_numParticles; ++i)

{

if (m_bFp64Positions)

glVertex3dv(&m_pos_fp64[k]);

else

{

glVertex3fv(&m_pos[k]);

}

k += 4;

}

}

glEnd();

}

else

{

glEnableClientState(GL_VERTEX_ARRAY);

glBindBufferARB(GL_ARRAY_BUFFER_ARB, m_pbo);

if (m_bFp64Positions)

glVertexPointer(4, GL_DOUBLE, 0, 0);

else

glVertexPointer(4, GL_FLOAT, 0, 0);

if (color)

{

glEnableClientState(GL_COLOR_ARRAY);

glBindBufferARB(GL_ARRAY_BUFFER_ARB, m_vboColor);

//glActiveTexture(GL_TEXTURE1);

//glTexCoordPointer(4, GL_FLOAT, 0, 0);

glColorPointer(4, GL_FLOAT, 0, 0);

}

glDrawArrays(GL_POINTS, 0, m_numParticles);

glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);

glDisableClientState(GL_VERTEX_ARRAY);

glDisableClientState(GL_COLOR_ARRAY);

}

}

void ParticleRenderer::display(DisplayMode mode /* = PARTICLE_POINTS */)

{

switch (mode)

{

case PARTICLE_POINTS:

glColor3f(1, 1, 1);

glPointSize(m_pointSize);

glUseProgram(m_programPoints);

_drawPoints();

glUseProgram(0);

break;

case PARTICLE_SPRITES:

default:

{

// punct de setare

glEnable(GL_POINT_SPRITE_ARB);

glTexEnvi(GL_POINT_SPRITE_ARB, GL_COORD_REPLACE_ARB, GL_TRUE);

glEnable(GL_VERTEX_PROGRAM_POINT_SIZE_NV);

glPointSize(m_spriteSize);

glBlendFunc(GL_SRC_ALPHA, GL_ONE);

glEnable(GL_BLEND);

glDepthMask(GL_FALSE);

glUseProgram(m_programSprites);

GLuint texLoc = glGetUniformLocation(m_programSprites, "splatTexture");

glUniform1i(texLoc, 0);

glActiveTextureARB(GL_TEXTURE0_ARB);

glBindTexture(GL_TEXTURE_2D, m_texture);

glColor3f(1, 1, 1);

glSecondaryColor3fv(m_baseColor);

_drawPoints();

glUseProgram(0);

glDisable(GL_POINT_SPRITE_ARB);

glDisable(GL_BLEND);

glDepthMask(GL_TRUE);

}

break;

case PARTICLE_SPRITES_COLOR:

{

// punct de setare

glEnable(GL_POINT_SPRITE_ARB);

glTexEnvi(GL_POINT_SPRITE_ARB, GL_COORD_REPLACE_ARB, GL_TRUE);

glEnable(GL_VERTEX_PROGRAM_POINT_SIZE_NV);

glPointSize(m_spriteSize);

glBlendFunc(GL_SRC_ALPHA, GL_ONE);

glEnable(GL_BLEND);

glDepthMask(GL_FALSE);

glUseProgram(m_programSprites);

GLuint texLoc = glGetUniformLocation(m_programSprites, "splatTexture");

glUniform1i(texLoc, 0);

glActiveTextureARB(GL_TEXTURE0_ARB);

glBindTexture(GL_TEXTURE_2D, m_texture);

glColor3f(1, 1, 1);

glSecondaryColor3fv(m_baseColor);

_drawPoints(true);

glUseProgram(0);

glDisable(GL_POINT_SPRITE_ARB);

glDisable(GL_BLEND);

glDepthMask(GL_TRUE);

}

break;

}

SDK_CHECK_ERROR_GL();

}

const char vertexShaderPoints[] =

{

"void main() \n"

"{ \n"

" vec4 vert = vec4(gl_Vertex.xyz, 1.0); \n"

" gl_Position = gl_ProjectionMatrix * gl_ModelViewMatrix * vert; \n"

" gl_FrontColor = gl_Color; \n"

"} \n"

};

const char vertexShader[] =

{

"void main() \n"

"{ \n"

" float pointSize = 500.0 * gl_Point.size; \n"

" vec4 vert = gl_Vertex; \n"

" vert.w = 1.0; \n"

" vec3 pos_eye = vec3 (gl_ModelViewMatrix * vert); \n"

" gl_PointSize = max(1.0, pointSize / (1.0 – pos_eye.z)); \n"

" gl_TexCoord[0] = gl_MultiTexCoord0; \n"

//" gl_TexCoord[1] = gl_MultiTexCoord1; \n"

" gl_Position = gl_ProjectionMatrix * gl_ModelViewMatrix * vert; \n"

" gl_FrontColor = gl_Color; \n"

" gl_FrontSecondaryColor = gl_SecondaryColor; \n"

"} \n"

};

const char pixelShader[] =

{

"uniform sampler2D splatTexture; \n"

"void main() \n"

"{ \n"

" vec4 color2 = gl_SecondaryColor; \n"

" vec4 color = (0.6 + 0.4 * gl_Color) * texture2D(splatTexture, gl_TexCoord[0].st); \n"

" gl_FragColor = \n"

" color * color2;\n"//mix(vec4(0.1, 0.0, 0.0, color.w), color2, color.w);\n"

"} \n"

};

void ParticleRenderer::_initGL()

{

m_vertexShader = glCreateShader(GL_VERTEX_SHADER);

m_vertexShaderPoints = glCreateShader(GL_VERTEX_SHADER);

m_pixelShader = glCreateShader(GL_FRAGMENT_SHADER);

const char *v = vertexShader;

const char *p = pixelShader;

glShaderSource(m_vertexShader, 1, &v, 0);

glShaderSource(m_pixelShader, 1, &p, 0);

const char *vp = vertexShaderPoints;

glShaderSource(m_vertexShaderPoints, 1, &vp, 0);

glCompileShader(m_vertexShader);

glCompileShader(m_vertexShaderPoints);

glCompileShader(m_pixelShader);

m_programSprites = glCreateProgram();

glAttachShader(m_programSprites, m_vertexShader);

glAttachShader(m_programSprites, m_pixelShader);

glLinkProgram(m_programSprites);

m_programPoints = glCreateProgram();

glAttachShader(m_programPoints, m_vertexShaderPoints);

glLinkProgram(m_programPoints);

_createTexture(32);

glGenBuffers(1, (GLuint *)&m_vboColor);

glBindBuffer(GL_ARRAY_BUFFER_ARB, m_vboColor);

glBufferData(GL_ARRAY_BUFFER_ARB, m_numParticles * 4 * sizeof(float), 0, GL_STATIC_DRAW_ARB);

glBindBuffer(GL_ARRAY_BUFFER_ARB, 0);

}

//Evaluare functii – denumire EvalHermite

inline float evalHermite(float pA, float pB, float vA, float vB, float u)

{

float u2=(u*u), u3=u2*u;

float B0 = 2*u3 – 3*u2 + 1;

float B1 = -2*u3 + 3*u2;

float B2 = u3 – 2*u2 + u;

float B3 = u3 – u;

return (B0*pA + B1*pB + B2*vA + B3*vB);

}

unsigned char *createGaussianMap(int N)

{

float *M = new float[2*N*N];

unsigned char *B = new unsigned char[4*N*N];

float X,Y,Y2,Dist;

float Incr = 2.0f/N;

int i=0;

int j = 0;

Y = -1.0f;

//float mmax = 0;

for (int y=0; y<N; y++, Y+=Incr)

{

Y2=Y*Y;

X = -1.0f;

for (int x=0; x<N; x++, X+=Incr, i+=2, j+=4)

{

Dist = (float)sqrtf(X*X+Y2);

if (Dist>1) Dist=1;

M[i+1] = M[i] = evalHermite(1.0f,0,0,0,Dist);

B[j+3] = B[j+2] = B[j+1] = B[j] = (unsigned char)(M[i] * 255);

}

}

delete [] M;

return (B);

}

void ParticleRenderer::_createTexture(int resolution)

{

unsigned char *data = createGaussianMap(resolution);

glGenTextures(1, (GLuint *)&m_texture);

glBindTexture(GL_TEXTURE_2D, m_texture);

glTexParameteri(GL_TEXTURE_2D, GL_GENERATE_MIPMAP_SGIS, GL_TRUE);

glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);

glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, resolution, resolution, 0,

GL_RGBA, GL_UNSIGNED_BYTE, data);

}

ANEXA 2. Lista figurilor si a tabelelor

Lista figuri :

Figura 1.1: Modul firelor de execuție.

Figura 2.1. Nuclee CPU – GPU

Figura 2.2. Calcul înmulțire matrice pe CPU

Figura 2.3. Calcul înmulțire matrice pe GPU – exemplificat

Figura 2.4. Reprezentarea timpului de execuție a programului “Inmultirea matricilor” folosind CPU, respectiv GPU

Figura 2.5. Reprezentarea timpului de execuție a programului “Permutari CUDA C” folosind CPU, respectiv GPU

Figura 3.1. Fereastră de la o randare interactivă 3D a unui sistem de 10240 corpuri, simulate prin aplicație

Figura 3.2. Figura schematica a unei placi de calcul

Figura 3.3. Exemplificare bloc fire de executie – placi calcul

Figura 3.4. Grila de blocuri care calculeaza toate fortele

Figura 3.5. Imagine de la rularea programului – prezentare performante

Figura 3.6. Cresterea performantei folosind multiple fire de executie per corp

Figura 3.7. Fereastra de la o randare interactiva 3D a unui sistem de 10240 corpuri, simulate prin aplicatie folosind precizie dubla

Lista tabele:

Tabelul 1.1. Citirea și scrierea – Memorie

Tabelul 2.1. Date acumulate pentru problema Înmulțirea Matricilor

Tabelul 2.2. Date acumulate pentru problema Permutări CUDA C

BIBLIOGRAFIE

Jason Sanders and Edward Kandrot – An introduction to general-purpose GPU programming

Shane Cook – “Cuda programming – A developer’s guide to parallel computing with GPU’s”

Aarseth S. – Gravitational N-Body Simulations

NVIDIA Corporation 2007 – Nvidia CUDA Compute Unified Device Arhitecture Programming Guide

Stephen Jones – Introduction to Dynamic Parallelism

Cuda Tutorial – https://developer.nvidia.com/cuda-zone

Prof. Mike Giles – CPUs and GPUs – http://people.maths.ox.ac.uk/gilesm/cuda/lecs/lec0.pdf

Parallel computing – https://en.wikipedia.org/?title=Parallel_computing

Blaise Barney – Introduction to Parallel Computing

Lars Nyland, Mark Harris and Jan Prins – Fast N-Body Simulation with CUDA

Sverre J. Aarseth – Gravitational N-Body Simulations

Introduction to CUDA C – http://www.nvidia.com/content/GTC-2010/pdfs/2131_GTC2010.pdf