Show page

Differences

This shows you the differences between two versions of the page.

--- asc:laboratoare:09 [2020/04/10 22:17]
florin.mihalache [Exercitii]
+++ asc:laboratoare:09 [2024/02/29 13:09] (current)
giorgiana.vlasceanu created
@@ Line 1: / Line 1: @@
-====== Laboratorul 09 - Advanced CUDA ======
+====== Laboratorul 09 - Tehnici de Optimizare de Cod – Inmultirea Matricelor ======
-===== Spatiu unificat memorie =====
+===== Obiective =====
-De la [[http://developer.download.nvidia.com/compute/cuda/6_0/rel/docs/CUDA_Toolkit_Release_Notes.pdf| CUDA 6.0]], NVIDIA a schimbat semnificativ modelul de programare prin facilitarea comunicarii unitatii CPU (host) cu unitatea GPU (device), in mod transparent prin acelasi set de adrese de memorie virtuale. Astfel exista posibilitatea ca prin acelasi pointer de memorie sa se scrie date atat de catre CPU cat si de catre GPU. Evident transferurile de memorie au loc intre spatii diferite de adresare (ex RAM vs VRAM), dar acest lucru se intampla transparent la nivel de aplicatie CUDA / pentru programator.
+In acest laborator vom exemplifica o serie de optimizari de cod pe una dintre cele mai simple, si in acelasi timp utilizate probleme, si anume, inmultirea matricelor.
-{{:asc:lab9:nv-unified.png?640|NVIDIA Unified Memory}}
+==== De ce inmultirea matricelor? ====
-Mai jos avem un exemplu de folosire a memoriei unificate. Singura diferenta fata de alocarea pe CPU/HOST este ca memoria trebuie alocata cu cudaMallocManaged si dealocata cu cudaFree.
+Este o operatie fundamentala si elementara in algebra liniara ce serveste la rezolvarea unui numar extrem de mare de probleme, cum ar fi: rezolvarea sistemelor liniare de ecuatii in majoritatea domeniilor stiintifice si economice (operatiile cu matrice sunt practic prezente pretudindeni); calcule si operatii cu grafuri; inversari de matrice. Problema inmultirii matricelor este in mod cert cea mai bine studiata problema in HPC (High Performance Computing), ea beneficiind de o multitudine de algoritmi inteligenti si implementari performante pe toate arhitecturile existente astazi. Pentru a simplifica lucrurile, in acest laborator ne vom ocupa doar de inmultirea matricelor patratice.
-<code C>
+===== Cel mai simplu algoritm =====
-#include <iostream>
-#include <math.h>
+Intuitiv, cel mai simplu algoritm, urmeaza formularea matematica:
+{{:asc:lab5:cij.jpg|}}
-// CUDA kernel to add elements of two arrays
+Matricele A = [aij], i,j=1,...,N  si B = [bij], i,j=1,...,N sunt salvate ca vectori bidimensionali de marime N x N. Matricea rezultat C = A x B = [cij], i,j=1,...,N, avand fireste aceeasi dimensiune.
-__global__
-void add(int n, float *x, float *y)
+{{:asc:lab5:axb_c.jpg|}}
-{
-  int index = blockIdx.x * blockDim.x + threadIdx.x;
+Cum este si de asteptat, similar cu majoritatea operatiilor din algebra liniara, formula de mai sus se transforma in urmatorul program extrem de simplu:
-  int stride = blockDim.x * gridDim.x;
-  for (int i = index; i < n; i += stride)
+<code cpp>
-    y[i] = x[i] + y[i];
+int i,j,k;
-}
+double a[N][N], b[N][N], c[N][N];
+// initializarea matricelor a si b
-int main(void)
+for (i=0;i<N;i++){
-{
+   for (j=0;j<N;j++){
-  int N = 1<<20;
+      c[i][j] = 0.0;
-  float *x, *y;
+      for (k=0;k<N;k++){
+	c[i][j] += a[i][k] * b[k][j];
-  // Allocate Unified Memory -- accessible from CPU or GPU
+      }
-  cudaMallocManaged(&x, N*sizeof(float));
+   }
-  cudaMallocManaged(&y, N*sizeof(float));
-  // initialize x and y arrays on the host
-  for (int i = 0; i < N; i++) {
-    x[i] = 1.0f;
-    y[i] = 2.0f;
-  }
-  // Launch kernel on 1M elements on the GPU
-  int blockSize = 256;
-  int numBlocks = (N + blockSize - 1) / blockSize;
-  add<<<numBlocks, blockSize>>>(N, x, y);
-  // Wait for GPU to finish before accessing on host
-  cudaDeviceSynchronize();
-  // Check for errors (all values should be 3.0f)
-  float maxError = 0.0f;
-  for (int i = 0; i < N; i++)
-    maxError = fmax(maxError, fabs(y[i]-3.0f));
-  std::cout << "Max error: " << maxError << std::endl;
-  // Free memory
-  cudaFree(x);
-  cudaFree(y);
-  return 0;
 }
 </code>
-===== Operatii atomice CUDA =====
+<note important>
-CUDA ofera acces la multiple operatii atomice tip citire-modificare-scriere. Acestea presupun serializarea accesului in contextul mai multor thread-uri. Functiile sunt limitate la anumite tipuri de date:
+**Cat de bun este acest algoritm?**
-  - int
-  - unsigned int
-  - unsigned long long int
-  - float
-  - double
-Exemple de functii atomice:
+Algoritmul este bun pentru ca:
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicadd|atomicAdd]]
+  *	Se poate specifica in doar cateva linii;
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicsub|atomicSub]]
+  *	Este o mapare directa a formulei de calcul pentru Cij (din algebra liniara); este usor de inteles si de urmarit de catre oricine poseda cunostinte minime de matematica;
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicexch|atomicExch]]
+  *	In sfarsit, in mod sigur nu contine bug-uri datorita simplitatii extreme pe care o manifesta algoritmul!
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicmin|atomicMin]]
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicmax|atomicMax]]
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicinc|atomicInc]]
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicdec|atomicDec]]
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicand|atomicAnd]]
-  - [[https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#atomicor|atomicOr]]
-<note important>
+Algoritmul este prost pentru ca:
-A se consulta cu atentie documentatia CUDA inainte de folosirea unei operatii atomice (legat de contextul in care se aplica, cum opereaza, limitari etc).
+  *	Are performante extrem de reduse!
-</note>
-In codul de mai jos se lanseaza un kernel concurrentRW in configuratie numBlocks=8 fiecare cu cate 10 thread-uri.
+De aceea ne vom ocupa in acest laborator de optimizarea acestei operatii din punctul de vedere al performantei.
-<code C>
+</note>
-#include <iostream>
-#define NUM_ELEM        8
+===== Optimizarea algoritmului de inmultire a doua matrice =====
-#define NUM_THREADS     10
-using namespace std;
+==== Detectarea constantelor din bucle ====
-__global__ void concurrentRW(int *data) {
+Prima optimizare, consta in a observa ca c[i][j] este o constanta in cadrul ciclului interior k. Totusi, pentru un compilator acest fapt nu este neaparat evident deoarece c[i][j] este o referinta in cadrul unui vector. Astfel, o prima optimizare va arata asa:
-...
+<code cpp>
+for (i=0;i<N;i++){
+   for (j=0;j<N;j++){
+      register double suma = 0.0;
+      for (k=0;k<N;k++) {
+         suma += a[i][k] * b[k][j];
+      }
+      c[i][j] = suma;
+   }
 }
+</code>
-int main(int argc, char *argv[]) {
+In acest mod, compilatorul va putea avea grija ca variabila suma sa fie tinut intr-un registru, permitand astfel o utilizare optima a acestei resurse. Astfel utilizarea keyword-ului "register" este util de folosit ca hint pentru compilator, atunci cand socotiti ca acest lucru este util.
-    int* data = NULL;
-    bool errorsDetected = false;
-    cudaMallocManaged(&data, NUM_ELEM * sizeof(unsigned long long int));
+==== Accesul la vectori ====
-    if (data == 0) {
-        cout << "[HOST] Couldn't allocate memory\n";
-        return 1;
-    }
-    // init all elements to 0
+Un alt aspect care necesita resurse din plin, este utilizarea si accesul variabilelor de tip vectorial. De fiecare data cand programul face o referinta la un obiect de tipul X[i][j][k] compilatorul trebuie sa genereze expresii aritmetice complexe, pentru a calcula aceasta adresa, in cadrul vectorului muldimensional X. De exemplu, iata cum arata un vector bidimensional in limbajul C (salvat row-major):
-    cudaMemset(data, 0, NUM_ELEM);
-    // launch kernel writes
+{{:asc:lab5:aij.jpg|}}
-    concurrentRW<<<NUM_ELEM, NUM_THREADS>>>(data);
-    cudaDeviceSynchronize();
-    if (cudaSuccess != cudaGetLastError()) {
-        return 1;
-    }
-    for(int i = 0; i < NUM_ELEM; i++) {
+Astfel, pentru ''N = 6, M = 4: a[2][3] = a[0][0] + 2*6 + 3 = a[0][0] + 15''
-        cout << i << ". " << data[i] << endl;
-        if(data[i] != (NUM_THREADS * (NUM_THREADS - 1) / 2)) {
-            errorsDetected = true;
-        }
-    }
-    if(errorsDetected) {
+In limbaje de programare ca FORTRAN-ul, formula este inversata, deoarece aceste limbaje salvează vectorii în format column-major:
-        cout << "Errors detected" << endl;
-    } else {
-        cout << "OK" << endl;
-    }
-    return 0;
+''a[i][j] = a[0][0] + j*M + i''
-}
-</code>
-Functie concurrentRW citeste valoarea de la adresa data[blockIdx.x], o incrementeaza cu threadId si apoi o scrie.
+Oricare ar fi asezarea vectorilor in memorie, accesele la vectori sunt scumpe din punctul de vedere al performantelor. Noi vom considera de aici inainte o asezare row-major, ca in limbajul C. Conform acestei formule, pentru vectori bidimensionali (matrice), fiecare acces presupune doua adunari si o inmultire (de numere intregi). Evident, pentru vectori cu mai multe dimensiuni, aceste costuri cresc considerabil. Astfel, in momentul in care compilatorul intalneste instructiunea:
-In acest caz avem 10 thread-uri care fac operatii citire/scriere la aceeasi adresa, deci un comportament nedefinit.
-<code C>
+''suma += a[i][k] * b[k][j]''
-__global__ void concurrentRW(int *data) {
-    // NUM_THREADS try to read and write at same location
+se vor efectua implicit, suplimentar inmultirii si adunarii in virgula mobila implicata de codul de mai sus, patru adunari si doua inmultiri in numere intregi pentru a calcula adresele necesare din vectorii a si b. Se intampla astfel destul de frecvent ca procesorul sa nu aiba date disponibile pentru a lucra in continuu, din cauza faptului ca overhead-ul pentru calculul adreselor este semnificativ.
-    data[blockIdx.x] = data[blockIdx.x] + threadIdx.x;
-}
+Astfel, un mod de a spori viteza programului este renuntarea la accesele vectoriale prin derefentiere utilizand in acest scop pointeri. De exemplu:
-</code>
-Exemplu rezultat:
+<code cpp>
-<code sh>
+for (j=0;j<N;j++)
-. 9
+   a[i][j] = 2;         // 2*N adunari si N inmultiri
-. 9
-. 9
-. 9
-. 9
-. 9
-. 9
-. 9
-Errors detected
 </code>
-Corect ar fi folosirea functie atomicAdd pentru a serializa accesul.
+se va inlocui cu:
-<code C>
+<code cpp>
-__global__ void concurrentRW(int *data) {
+double *ptr=&(a[i][0]); // 2 adunari si o inmultire
-    // NUM_THREADS try to read and write at same location
+for (j=0;j<N;j++) {
-    atomicAdd(&data[blockIdx.x], threadIdx.x);
+   *ptr = 2;
+   ptr++;               // N adunari in numere intregi
 }
 </code>
-Rezultatul rularii este:
+In mod similar se procedeaza si pentru cazul in care indexul incrementat este cel al liniilor si nu cel al coloanelor. In ambele cazuri, practic se va calcula "de mana" adresa in cadrul vectorului, exact in modul in care ar face-o compilatorul limbajului folosit. Totusi, rezolvarea noastra este mai rapida, deoarece ea tine cont de pozitia in care ne aflam in cadrul vectorului, lucru destul de complicat de facut automat. De exemplu, pentru a trece la urmatoarea coloana, e suficient sa adunam N pointer-ului, fata de recalcularea pornind de la @(a[0][0]) ce necesita doua inmultiri si o adunare in intregi. Evident, facilitatile oferite de limbaje ca C-ul, ne vin in ajutor: astfel incrementarile de pointeri de tip char * vor face incrementarea cu un byte, in vreme ce pentru int * se va face cu patru bytes. Ca urmare a aspectelor prezentate mai sus, iata forma optimizata in care ajunge algoritmul nostru:
-<code sh>
-. 45
+<code cpp>
-. 45
+for(i = 0; i < N; i++){
-. 45
+  double *orig_pa = &a[i][0];
-. 45
+  for(j = 0; j < N; j++){
-. 45
+    double *pa = orig_pa;
-. 45
+    double *pb = &b[0][j];
-. 45
+    register double suma = 0;
-. 45
+    for(k = 0; k < N; k++){
-OK
+      suma += *pa * *pb;
+      pa++;
+      pb += N;
+    }
+    c[i][j] = suma;
+  }
+}
 </code>
-==== Operatii atomice system wide ====
+<note tip>Atentie! Codul de mai sus va da rezultate corecte doar daca matricile sunt declarate global sau pe stivă pentru că în felul acesta sunt stocate continuu în memorie (și are sens pb += N). Dacă alocați dinamic, atunci folosiți matrici liniarizate și adaptați acest cod pentru cazul lor.</note>
-Unitatile GPU ce au Compute capability 6.x permit largirea scopului operatiilor atomice. De exemplu atomicAdd_system garanteaza ca operatia este atomica cand atat thread-urile de pe unitatea GPU cat si cele de pe unitatea CPU incearca sa acceseze datele. Mai jos avem un exemplu de folosire al functiei atomicAdd_system.
+<note important>
+Din primele doua optimizari se pot desprinde cateva concluzii. Prima ar fi ca optimizarea unui cod (din punct de vedere al performantelor), presupune utilizarea a cat mai putine constructii complexe (high-level), puse la dispozitie de limbajul folosit. Aceasta concluzie poate suna extrem de ciudat pentru cineva care porneste de la ideea ca facilitatile limbajelor de programare sunt acolo pentru a fi folosite. Da, este adevarat acest lucru, insa atunci cand vrei performanta, trebuie sa stii ce constructii sa eviti! Astfel, apare concluzia a doua: vectorii sunt concepte mai abstracte decat pointerii (ca implementare), asadar, utilizati pointeri cand vreti viteza. Viteza crescuta insa, va fi obtinuta cu pretul unui cod mult mai dificil de urmarit si de inteles, mai rau, mult mai greu de debug-at. Un cod complex si performant, de multe ori poate contine bug-uri extrem de subtile si greu de depistat. Asadar, e util sa stii exact ceea ce faci cand incepi sa faci astfel de optimizari!
-<code C>
+</note>
-__global__ void mykernel(int *addr) {
-  atomicAdd_system(addr, 10);       // only available on devices with compute capability 6.x
+----
+==== Activitate practica - Optimizare constantelor si al accesului la vectori ====
+Intrebarea este acum: aduc ceva imbunatatiri optimizarile 1 si 2? Pentru a afla raspunsul la aceasta intrebare, va invitam sa implementati problema, cu optimizarile sugerate, si sa observati singuri ce se intampla.
+----
+==== Optimizarea pentru accesul la memorie ====
+Dupa cum ar trebui sa va fie destul de evident pana acum, din experienta voastra de programatori, memoria este in general cel mai problematic bottleneck. Optimizarile prezentate mai sus reduc timpul de executie intr-o oarecare masura, insa ele nu schimba in nici un fel modul in care memoria este accesata in cadrul algoritmului. Cu alte cuvinte, aceleasi locatii de memorie sunt accesate in aceeasi ordine, indiferent daca am operat sau nu optimizarile prezentate. O intrebare interesanta ar fi acum: ce se intampla, daca am schimba ordinea in care se executa buclele? S-ar obtine performante diferite?
+Pentru problema noastra, care contine trei bucle, exista asadar sase secvente posibile, si anume: i-j-k, i-k-j, j-i-k, j-k-i, k-i-j, si k-j-i. Fiecare dintre aceste secvente corespunde unui tip diferit de acces la memorie pentru matricele considerate. Deoarece bucla interioara este cea mai des executata, ne vom concentra acum atentia un pic asupra ei. Operatia executata acolo ramane:
+''c[i][j] += a[i][k] * b[k][j]''
+Pentru fiecare dintre cele trei matrice, a, b si c, fiecare element poate fi accesat in trei moduri diferite, si anume:
+  * Constant: accesul nu depinde de indexul buclei interioare
+  * Secvential: accesul la memorie este contiguu (adica in celule succesive de memorie)
+  * Nesecvential: accesul la memorie nu este contiguu (celulele de memorie logic succesive, sunt de fapt adresate cu pauze de dimensiune N)
+Astfel, pentru cele sase configuratii, se obtine:
+^ Loop order ^ c[i][j] += ^ a[i][k] ^ * b[k][j] ^
+| i-j-k: | Constant   | Secvential | Nesecvential |
+| i-k-j: | Secvential | Constant | Secvential |
+| j-i-k: | Constant | Secvential | Nesecvential |
+| j-k-i: | Nesecvential | Nesecvential | Constant |
+| k-i-j: | Secvential | Constant | Secvential |
+| k-j-i: | Nesecvential | Nesecvential | Constant |
+Care sunt totusi, comparativ, performantele celor trei moduri de acces? In mod clar, accesul constant este mai bun decat cel secvential – aceste constante in cadrul unor bucle, sunt in general puse in registri, ducand la imbunatatirea performantelor  algoritmului, dupa cum s-a aratat in optimizarea 1. Accesul secvential la randul sau, este mai bun decat cel nesecvential, in principal pentru ca utilizeaza considerabil mai bine cache-ul.
+<note important>
+Luand in considerare aceste observatii, putem concluziona ca:
+  * Configuratiile k-i-j si i-k-j ar trebui sa aiba cele mai bune performante
+  * Configuratiile i-j-k si j-i-k ar trebui sa fie mai proaste decat primele, si
+  * Configuratiile j-k-i si k-j-i ar trebui sa fie cele mai proaste!
+</note>
+==== Activitate practica - Ordinea buclelor ====
+----
+Efectiv, care este adevarul? Construiti singuri aceste scenarii si analizati aceasta problema!
+----
+Pentru a studia mai in detaliu problema, sa analizam un pic configuratia i-j-k (desi nu este cea mai buna configuratie, cum vedem de mai sus):
+<code cpp>
+for (i=0;i<N;i++){
+   for (j=0;j<N;j++){
+      sum=0;
+      for (k=0;k<N;k++)
+         sum+=a[i][k]*b[k][j];
+      c[i][j] = sum;
+   }
 }
+</code>
-void foo() {
+Cate cache-miss-uri sunt generate in acest algoritm, cu aceasta secventa de acces la memorie? In mod evident, aceasta nu este o intrebare usoara. De exemplu: daca fiecare matrice ar fi de doua ori mai mare decat cache-ul, ar avea loc multe incarcari si eliberari de linii, ducand astfel la o formula complicata. Astfel, cel mai simplu aproximam, si consideram ca dimensiunea matricei este mult mai mare decat cea a Cache-ului. Astfel, fie C, numarul de elemente din matrice ce intra in Cache.
-  int *addr;
+Astfel, considerand algoritmul de mai sus (fara optimizarea pentru constante):
-  cudaMallocManaged(&addr, 4);
-  *addr = 0;
-   mykernel<<<...>>>(addr);
+<code cpp>
-   __sync_fetch_and_add(addr, 10);  // CPU atomic operation
+for (i=0;i<N;i++){
+   // Citeste linia i pt a in Cache (Ra)
+   // Scrie linia i a lui c in Memorie (Wc)
+   for (j=0;j<N;j++){
+      // Citeste coloana j a lui b in Cache (Rb)
+      for (k=0;k<N;k++){
+         c[i][j] += a[i][k] * b[k][j];
+      }
+   }
 }
 </code>
-===== Operatii asincrone CUDA =====
+Astfel, daca L este dimensiunea unei linii de Cache: pentru (Ra) obtinem aproximativ N*(N/L) cache-miss-uri, pentru (Wc) la fel, iar pentru (Rb) un dezastruos N*N*N! Acest lucru se intampla deoarece, desi accesul la b este secvential intr-o coloana, matricea este salvata in memorie utilizand row-major! Concluzia este descurajatoare: 2N<sup>2</sup>/L + N<sup>3</sup> -> N<sup>3</sup> cache-miss-uri! Se adauga la acest aspect si cele 2N<sup>3</sup> operatii aritmetice, si se ajunge la raportul: operatii aritmetice / operatii cu memoria -> 2. Acest lucru este extrem de rau, deoarece noi stim de la (curs) si de la alte materii, ca arhitecturile calculatoarelor NU sunt echilibrate, si ca operatiile aritmetice sunt de ordine de marime mai rapide decat operatiile cu memoria. De aceea, memoria ramane in continuare bottleneck-ul pentru aceasta implementare a inmultirii de matrice. Pentru a obtine performante mai bune, este necesara obtinerea unui raport considerabil mai mare.
-In CUDA, urmatoarele operatii sunt definite ca fiind independente si pot fi executate concurent:
+Cum se face insa, ca pentru N<sup>2</sup> elemente intr-o matrice, ajungem la N<sup>3</sup> cache-miss-uri? Pai am stabilit ca acest lucru se datoreaza accesului ineficient al lui b, deoarece se incearca incarcarea coloana cu coloana a matricei!
-  - Calcule pe unitatea host
-  - Calcule pe unitatea device
-  - Transfer memorie host -> device
-  - Transfer memorie device -> host
-  - Transfer memorie device -> device
-Nivelul de concurenta o sa depinda si de capabilitatea unitatilor GPU (compute capability).
-In continuare vom explora mai multe scenarii de executie concurenta a operatiilor descrise.
-==== Executie asincrona Host si Device ====
+Concluzia acestei analize este ca nu putem spune, doar dupa numarul de operatii efectuate si dimensiunea datelor folosite, daca un algoritm va suferi sau nu din cauza unui bottleneck la memorie.
-Folosind apeluri asincrone, operatiile de executie catre device sunt puse in coada avand controlul intors catre host instant. Astfel unitatea host poate continua executia fara sa fie blocata in asteptarea executiei.
+Solutia este: utilizarea mai ingenioasa a cache-ului.
-Urmatoarele operatii sunt asincrone relativ la host:
-  - Lansari de kernel
-  - Copieri in cadrul spatiului de memorie a unui device
-  - Copiere memorie host -> device, avand < 64 KB
-  - Copiere memorie host -> device, avand functii cu sufix Async
-  - Functii memorie set
-Pentru a face debug unor scenarii de executie asincrona se poate dezactiva complet executia asincrona setand variabila de mediu CUDA_LAUNCH_BLOCKING la 1. Executia de kernels este sincrona cand se ruleaza cu un profiler (Nsight, Visual Profiler).
+Acest lucru se poate realiza prin reorganizarea operatiilor din cadrul inmultirii de matrice pentru a obtine mai multe cache-hit-uri. Faptul ca adunarea si inmultirea sunt atat operatii asociative, cat si comutative face posibila aceasta reordonare a operatiilor. Acesta este un subiect de cercetare asupra caruia si-au indreptat atentia numerosi cercetatori de-a lungul timpului, generand o multitudine de algoritmi si de teoreme matematice care sa ii sustina. In orice caz, daca vom considera r = raportul intre operatiile aritmetice si operatiile  la memorie (cu cache-miss-uri), este evident ca se doreste un r maxim, pentru a elimina bottleneck-ul de la memorie. S-a aratat ca orice reorganizare a acestui algoritm este limitata la r = O(sqrt(C)), unde C este dimensiunea Cache-ului (in numar de elemente ce intra in Cache). Acest lucru arata ca r nu scaleaza cu dimensiunea matricei N, indiferent de impartirea intuitiva a lui 2N<sup>3</sup> la N<sup>2</sup>...
-<code C>
+==== Solutia: “Blocked Matrix Multiplication” ====
-result = cudaMemcpyAsync(d_a, a, N, cudaMemcpyHostToDevice, stream1)
+Pentru a rezolva problema accesului in b pentru coloane intregi, se va trece la accesarea unui subset a unei coloane in b, sau a mai multor coloane la un moment dat. Pentru o mai buna intelegere, urmariti desenele de mai jos:
+{{:asc:lab5:c_axb_1.jpg|}}
+Ideea de baza este refolosirea cat mai buna a elementelor aflate in cache (pentru matricea b). Astfel odata cu calculul lui c[i][j], de ce nu am calcula si c[i][j+1], daca tot se afla in cache si coloana j+1. Acest lucru presupune insa reordonarea operatiilor astfel: calculeaza primii b termeni pentru c[i][j], calculeaza primii b termeni pentru c[i][j+1], calculeaza urmatorii b termeni pentru c[i][j], calculeaza urmatorii b termeni pentru c[i][j+1], etc.
+{{:asc:lab5:c_axb_2.jpg|}}
+In acest mod, de ce nu am calcula o intreaga sectiune de linie din c, folosind aceste reordonari de operatii?
+Ce s-ar intampla daca am incerca sa calculam o intreaga linie din c?
+{{:asc:lab5:c_axb_3.jpg|}}
+Ar insemna ca trebuie sa incarcam toate coloanele lui b in memorie (cache), lucru pe care am incercat sa il evitam aici! Astfel, se vor refolosi doar acele blocuri din b ce au fost deja incarcate. De aici nu ne mai ramane decat sa utilizam intreaga linie de cache din b, si obtinem ideea de baza a algoritmului “Blocked Matrix Multiplication”:
+{{:asc:lab5:c_axb_4.jpg|}}
+Operatiile trebuie reordonate astfel: calculeaza primii b termeni pentru c[i][j] din blocul C, calculeaza urmatorii b termeni pentru c[i][j] din blocul C, ..., calculeaza ultimii b termeni pentru c[i][j] din blocul C. Generalizand:
+{{:asc:lab5:bmm.jpg|}}
+Pentru a calcula blocul C<sub>22</sub> folosim formula:
+C<sub>22</sub> %%=%% A<sub>21</sub>B<sub>12</sub> + A<sub>22</sub>B<sub>22</sub> + A<sub>23</sub>B<sub>32</sub> + A<sub>24</sub>B<sub>42</sub>
+ce presupune patru inmultiri si patru adunari de matrice. Ideea este ca fiecare inmultire opereaza pe un block suficient de mic ca dimensiune astfel incat sa intre in Cache!
+Versiunea inmultirii de matrice utilizand metoda bloc si ordonarea i-j-k devine:
+<code cpp>
+for (i=0;i<N/b;i++){
+   for (j=0;j<N/b;j++){
+      for (k=0;k<N/b;k++){
+         C[i][j] += A[i][k]*B[k][j]
+      }
+   }
+}
 </code>
-==== Executie asincrona programe kernel ====
+unde:
+  * b este dimensiunea blocului (presupunem ca b divide N)
+  * C[i][j] este un bloc al matricei C pe linia i si coloana j
+  * "+=" inseamna adunare de matrice
+  * si "*" inseamna inmultire de matrice
-Arhitecturile cu compute capability 2.x sau mai nou, pot executa in paralel instante de kernel diferite. Aceste unitati de executie o sa aibe proprietate concurrentKernels setata la 1 (se face query la device properties inainte). Numarul maxim de lansari asincrone de kernele diferite este dependent de arhitectura (se verifica in functie de compute capability). Singura restrictie este ca programele kernel sa fie in acelasi context.
+Ce se intampla cu Cache-miss-urile acum?
+<code cpp>
+for (i=0;i<N/b;i++){
+   for (j=0;j<N/b;j++){
+      // Scrie blocul C[i][j] al lui c in Memorie (Wc)
+      for (k=0;k<N/b;k++){
+         // Citeste blocul A[i][k] pt a in Cache (Ra)
+         // Citeste blocul B[k][j] pt b in Cache (Rb)
+         C[i][j] += A[i][k] * B[k][j];
+      }
+   }
+}
+</code>
-==== Executie si transfer date asincron  ====
+Pentru (Wc) avem acum (N/b)*(N/b)*b*b Cache-miss-uri, in vreme ce pentru (Ra) si (Rb) avem (N/b)*(N/b)*(N/b)*b*b, astfel ducand la N<sup>2</sup> + 2N<sup>3</sup>/b -> 2N<sup>3</sup>/b Cache-miss-uri pentru intregul algoritm. Combinand acest calcul cu faptul ca avem 2N<sup>3</sup> operatii aritmetice, rezulta un raport r = 2N<sup>3</sup>/b / 2N<sup>3</sup> -> b. Dupa cum am stabilit, r trebuie sa fie maxim (mai mare oricum decat 2-ul obtinut in varianta anterioara). Daca mergem pana la cazul extrem, il vom face pe b = N, dar asta nu este viabil, pentru ca atunci suntem din nou la cazul fara blocuri, de la care tocmai venim...
-Anumite device-uri pot executa un transfer asincron memorie alaturi de o executie de kernel. Acest lucru este dependent de compute capability si se poate verifica in device property asyncEngineCount. Deasemenea se pot face transferuri de memorie intra-device simultan cu executia de kernel cand atat device property concurrentKernels si asyncEngineCount sunt 1.
+Astfel, acest algoritm functioneaza doar daca blocurile intra in Cache. Acest lucru inseamna ca trei blocuri diferite, de dimensiune b*b, trebuie sa intre in Cache, pentru toate cele trei matrice (a, b si c). Daca C este dimensiunea Cache-ului in elemente de matrice, atunci trebuie sa fie 3b<sup>2</sup> ≤ C sau b ≤ √(C / 3) . Astfel, in cel mai bun caz, r-ul trebuie sa fie si el  √(C / 3).
-===== Dynamic Paralellism  =====
+Putem astfel spune, pentru diverse procesoare, cunoscand rata de operatii aritmetice la cache-miss-uri r, care este dimensiunea necesara a Cache-ului, pentru a rula acest algoritm, astfel incat procesorul sa NU astepte niciodata memoria:
-Paralelismul dinamic consta in posibilitatea de a lansa programe kernel din thread-urile ce ruleaza pe device/GPU. In alte cuvinte, unitatea GPU poate sa isi  atribuie noi task-uri/thread-uri fara interventia unitatii host/CPU. Aceasta manifestare este utila in problemele unde maparea threaduri<->date nu este simpla/triviala. De exemplu, in situatia unde unele thread-uri ar avea prea putin de lucru, iar altele prea mult (imaginea de mai jos, simulare fluide) - o situatia debalansata computational.
+/* TODO: PLEASE UPDATE THIS!
+^ Procesor ^ Dimensiune Cache (KB) ^
+|Ultra 2i   | 14.8 |
+|Ultra 3    | 4.7 |
+|Pentium 3  | 0.9 |
+|Pentium 3M | 2.4 |
+|Power 3    | 1.8 |
+|Power 4    | 5.4 |
+|Itanium 1  | 31.1 |
+|Itanium 2  | 0.7 |
-{{:asc:lab9:dynamic-paralellism.png?560|Fluid simulation}}
+*/
-Cerintele pentru paralelism dinamic sunt CUDA 5.0 ca Toolkit si respectiv Compute Capability 3.5. O lista cu GPU-uri NVIDIA si Compute Capability se regaseste [[https://developer.nvidia.com/cuda-gpus|aici]]. Pana acum am lucrat pe coada hp-sl.q care are GPU-uri K40M (Compute Capability 3.5) si ibm-dp.q/ibm48-dp.q care are GPU-uri C2050 (Compute Capability 2.0). Astfel doar coada hp-sl.q suporta paralelism dinamic pe placile GPU K40M.
+==== Activitate practica - BMM & Optimizare pentru Cache ====
-<code C>
+De aceea incercati sa experimentati cele prezentate in acest laborator, in C. Pentru cei interesati, incercati completarea tabelului de mai sus cu dimensiunea Cache-ului pentru procesoarele voastre personale. Acest lucru presupune evident, si o documentare asupra caracteristicilor sistemului propriu (determinarea r-ului, a dimensiunii Cache-ului etc.).
+<code cpp>
 #include <stdio.h>
-__global__ void childKernel()
+#include <stdlib.h>
+#include <sys/time.h>
+void BMMultiply(int n, double** a, double** b, double** c)
 {
-   printf("Hello ");
+    int bi=0;
+    int bj=0;
+    int bk=0;
+    int i=0;
+    int j=0;
+    int k=0;
+    // TODO: set block dimension blockSize
+    int blockSize=100;
+    for(bi=0; bi<n; bi+=blockSize)
+        for(bj=0; bj<n; bj+=blockSize)
+            for(bk=0; bk<n; bk+=blockSize)
+                for(i=0; i<blockSize; i++)
+                    for(j=0; j<blockSize; j++)
+                        for(k=0; k<blockSize; k++)
+                            c[bi+i][bj+j] += a[bi+i][bk+k]*b[bk+k][bj+j];
 }
-__global__ void parentKernel()
+int main(void)
 {
-   // launch child
+    int n;
-   childKernel<<<1,1>>>();
+    double** A;
-   if (cudaSuccess != cudaGetLastError()) {
+    double** B;
-      return;
+    double** C;
-   }
+    int numreps = 10;
-   // wait for child to complete
+    int i=0;
-   if (cudaSuccess != cudaDeviceSynchronize()) {
+    int j=0;
-      return;
+    struct timeval tv1, tv2;
-   }
+    struct timezone tz;
-   printf("World!\n");
+    double elapsed;
-}
+    // TODO: set matrix dimension n
-int main(int argc, char *argv[])
+    n = 500;
-{
+    // allocate memory for the matrices
-   // launch parent
-   parentKernel<<<1,1>>>();
+    // TODO: allocate matrices A, B & C
-   if (cudaSuccess != cudaGetLastError()) {
+    ///////////////////// Matrix A //////////////////////////
-      return 1;
+    // TODO ...
-   }
-   // wait for parent to complete
+    ///////////////////// Matrix B //////////////////////////
-   if (cudaSuccess != cudaDeviceSynchronize()) {
+    // TODO ...
-      return 2;
-   }
+    ///////////////////// Matrix C //////////////////////////
-   return 0;
+    // TODO ...
+    // Initialize matrices A & B
+    for(i=0; i<n; i++)
+    {
+        for(j=0; j<n; j++)
+        {
+            A[i][j] = 1;
+            B[i][j] = 2;
+        }
+    }
+    //multiply matrices
+    printf("Multiply matrices %d times...\n", numreps);
+    for (i=0; i<numreps; i++)
+    {
+        gettimeofday(&tv1, &tz);
+        BMMultiply(n,A,B,C);
+        gettimeofday(&tv2, &tz);
+        elapsed += (double) (tv2.tv_sec-tv1.tv_sec) + (double) (tv2.tv_usec-tv1.tv_usec) * 1.e-6;
+    }
+    printf("Time = %lf \n",elapsed);
+    //deallocate memory for matrices A, B & C
+    // TODO ...
+    return 0;
 }
 </code>
+===== In loc de concluzie =====
+Intelegerea reala a comportamentului unei aplicatii (algoritm), din punctul de vedere al utilizarii cache-ului (si al performantelor in general), este o chestiune complexa, ce necesita multa rabdare si cunostinte diverse. Deseori, aproximatii utile pot fi folosite pentru a imbunatati unele aspecte ale implementarii curente. Utilizarea blocurilor este intalnita deseori in algoritmi si aplicatii ce necesita performante crescute.
 ===== Exercitii =====
-  - logati-va pe ''fep.grid.pub.ro'' folosind contul de pe ''cs.curs.pub.ro''
-  - executati comanda ''wget https://ocw.cs.pub.ro/courses/_media/asc/lab9/lab9_skl.tar.gz -O lab9_skl.tar.gz''
-  - dezarhivati folosind comanda ''tar -xzvf lab9_skl.tar.gz''
-  - executati comanda ''qlogin -q hp-sl.q'' pentru a intra pe o statie specializata in calcul folosind GPU-uri
-  - incarcati modului de Nvidia CUDA folosind comanda ''module load libraries/cuda''
-<note tip>
+**Task 0**  - Rulați ''task01'', ''task02'' si ''task03''ca exemple pentru [[#activitate_practica_-_optimizare_constantelor_si_al_accesului_la_vectori | Optimizarea constantelor si al accesului la vectori]] folosind matrici liniarizate.
-Debug aplicatii CUDA [[https://docs.nvidia.com/cuda/cuda-gdb/index.html#introduction|aici]]
-</note>
-<note tip>
+**Task 1**  -  Implemenati  [[#activitate_practica_-_ordinea_buclelor | Ordonarea buclelor]] folosind matrici liniarizate in ''task11.c'' si ''task12.c''. Rulati si observati diferentele.
-Profiling aplicatii CUDA via nvprof [[https://docs.nvidia.com/cuda/profiler-users-guide/index.html#nvprof-overview|aici]]
-</note>
-  * Deschideți fișierul task1_unified.cu și completati cu alocarile/trasferurile de memorie cerute. (3p)
+**Task 2**  - Implementati  [[#activitate_practica_-_bmm_optimizare_pentru_cache | Optimizari pentru Cache]] folosind matrici liniarizate.
-    * In functia compute_NoUnifiedMem se va aloca memorie cu cudaMalloc.
+    * ''task21.c'' Implementati BMM
-    * In functia compute_UnifiedMem se va aloca memorie cu cudaMallocManaged.
+    * ''task22.c'' Implementati BMM alaturi de  [[#activitate_practica_-_optimizare_constantelor_si_al_accesului_la_vectori | Optimizarea constantelor si al accesului la vectori]].
-    * Folosind nvprof analizati cum ruleaza cele 2 programe:
-      * nvprof ./task1_mem_plain
-      * nvprof ./task1_mem_unified
-  * Deschideti fisierul task2_atomic.cu si urmăriți instrucțiunile TODO (3p)
-    * Completati functiile kernel_no_atomics, kernel_partial_atomics si kernel_full_atomics.
-    * Folosind nvprof analizati cum ruleaza cele 3 programe:
-      * nvprof ./task2_no_atomic
-      * nvprof ./task2_partial_atomic
-      * nvprof ./task2_full_atomic
-  * Deschideti fisierul task3_atomic_exch.cu si urmăriți instrucțiunile TODO (2p)
-    * Se da un vector cu sloturi ocupate (.raw!=0) sau libere (.raw=0) si un set de elemente de inserat.
-  * Deschideti fisierul task4_dynp.cu si urmăriți instrucțiunile TODO (2p)
-    * Folosind dynamic parallelism se va calcula suma primelor data[i] elemente din vectorul data[]
-<note important>
-Recomandăm sa va delogati mereu de pe serverele din cluster dupa terminarea sesiunii, utilizand comanda ''logout''
-Alternativ, daca ati uitat sesiuni deschise, puteti verifica acest lucru de pe fep.grid.pub.ro, utilizand comanda ''qstat''. In cazul in care identificati astfel de sesiuni "agatate", le puteti sterge (si va rugam sa faceti asta), utilizand comanda ''qdel -f ID'' unde ID-ul il identificati din comanda anterioara ''qstat''.
+**Task 3**  - In fisierul ''task3.c'' implementati [[#activitate_practica_-_ordinea_buclelor | Ordonarea buclelor]] i-k-j folosind [[#activitate_practica_-_optimizare_constantelor_si_al_accesului_la_vectori | Optimizarea constantelor si al accesului la vectori]]
-Daca nu veti face aceasta delogare, veti putea ajunge in situatia in care sa nu va mai puteti loga pe nodurile din cluster.
+**Task 4**  -  (Bonus) In general, nu recomandam alocarea matricelor ca vectori de vectori. Ca bonus va sugeram sa realizati un test unde se face acest tip de alocare si se verifica "performantele" obtinute.
+<note important>
+Toate rularile din acest laborator trebuiesc facute cu optiunea -O0 (O zero), adica ''gcc -O0 -o binary source-file.c''
+Motivul principal este ca optimizarile de compilator ar putea sa ascunda optimizarile pe care le veti incerca (si reusi) voi la laborator.
 </note>
 ===== Resurse =====
+  * Responsabilul acestui laborator: [[emil.slusanschi@cs.pub.ro|Emil Slușanschi]]
+  * <html><a class="media mediafile mf_pdf" href=":asc:lab5:index?do=export_pdf">PDF laborator</a></html>
+  * {{:asc:lab5:lab5_skl.tar.gz|Schelet Laborator 5}}
+  <hidden>  * {{:asc:lab5:sol:lab5_sol.tar.gz|Soluție Laborator 5}}   </hidden>
-{{:asc:lab9:lab9_skl.tar.gz|Schelet Laborator 9}}
-{{:asc:lab9:sol:lab9_sol.tar.gz|Solutie Laborator 9}}
+==== Discutii interesante ====
-{{:asc:lab9:asc_lab9.pdf|Enunt Laborator 9}}
+  * [[http://stackoverflow.com/questions/11227809/why-is-processing-a-sorted-array-faster-than-an-unsorted-array | De ce este mai rapida procesarea unui vector ordonat?]]
-  * Responsabili laborator: Grigore Lupescu, Vlad Ştefănescu, Tudor Barbu, Mihai Despotovici
+  * {{:asc:lab5:what_every_programmer_should_know_about_memory_by_ulrich_drepper_.pdf|What every programmer should know about memory.pdf}}
-==== Referinte ====
+==== Valgrind ====
+  * http://valgrind.org/docs/manual/cg-manual.html
-  * Documentatie CUDA:
+==== Referințe ====
-    * [[https://docs.nvidia.com/pdf/CUDA_C_Programming_Guide.pdf|CUDA C Programming]]
+<hidden>
-    * [[https://docs.nvidia.com/cuda/pdf/CUDA_Compiler_Driver_NVCC.pdf| CUDA NVCC compiler]]
+  * [[https://github.com/metallurgix | Exemple inmultire matrice]] <-- outdated
-    * [[https://docs.nvidia.com/cuda/profiler-users-guide/index.html| CUDA Visual Profiler]]
+</hidden>
-    * [[https://developer.download.nvidia.com/compute/cuda/9.1/Prod/docs/sidebar/CUDA_Toolkit_Release_Notes.pdf|CUDA 9.1 Toolkit]]
-    * [[https://developer.nvidia.com/cuda-gpus|CUDA GPUs]]
-  * Acceleratoare hp-sl.q (hpsl-wn01, hpsl-wn02, hpsl-wn03)
-    * [[http://international.download.nvidia.com/tesla/pdf/tesla-k40-passive-board-spec.pdf|NVIDIA Tesla K40M]]
-    * [[https://en.wikipedia.org/wiki/Nvidia_Tesla|NVIDIA Tesla]]
-  * Acceleratoare ibm-dp.q (dp-wn01, dp-wn02, dp-wn03)
-    * [[https://www.nvidia.com/docs/IO/43395/NV_DS_Tesla_C2050_C2070_jul10_lores.pdf|NVIDIA Tesla C2070]]
-    * [[http://www.nvidia.com/docs/io/43395/nv_ds_tesla_c2050_c2070_apr10_final_lores.pdf|NVIDIA Tesla 2050/2070]]
-    * [[https://cseweb.ucsd.edu/classes/fa12/cse141/pdf/09/GPU_Gahagan_FA12.pdf|NVIDIA CUDA Fermi/Tesla]]
-  * Advanced CUDA
-    * [[https://devblogs.nvidia.com/gpu-pro-tip-cuda-7-streams-simplify-concurrency/|CUDA Streams]]
-    * [[https://devblogs.nvidia.com/introduction-cuda-dynamic-parallelism/|CUDA Dynamic Parallelism]]