Presentazione a cura di:

1. Parallel sparse Matrix-Vectorand Matrix-Transpose-Vectormultiplication using compressed sparse blocks Presentazione a cura di: Marco Cherubini, Andrea De Pirro, David Santucci,Andrea Tersigni, Luca Tracuzzi A. Buluc, J. T. Fineman, M. Frigo, J. R. Gilbert, C. E. Leiserson Calcolo Parallelo e Distribuito Anno Accademico 2009/2010

2. Sommario • Formati di memorizzazione convenzionali • Il nuovo formato CSB • Moltiplicazione Matrice-Vettore con CSB • Analisi della complessità • Sperimentazione

3. Formati convenzionali: CSR • Analizziamo alcuni formati di memorizzazione convenzionali • Consideriamo matrici sparse n×n con nnz elementi non nulli • CSR - Compressed Sparse Rows • Memorizzazione per righe • Efficiente: memorizza n+nnz indici o puntatori • per matrici sparse n×n con nnz elementi non nulli • Adatto per y  Ax • Non adatto per y  ATx

4. Formati convenzionali: CSR

5. Ax parallelo con CSR   CSR_SPMV (A, x, y) 1     nA.rows 2     for i0 to n−1 inparallel 3         do y[i]0 4             for kA.row_ptr[i] to A.row_ptr[i+1]−1 5                 do y[i]y[i]+A.val[k]·x[A.col_ind[k]] • val[nnz] : array dei valori non nulli della matrice (ordinati per righe) • col_ind[nnz] : indici di colonna degli elementi nell'array val • row_ptr[n] : puntatori al'inizio della riga n nell'array val • Nota: ATx con CSC è analogo

6. Formati convenzionali: CSC • CSC - Compressed Sparse Columns • Memorizzazione per colonne • Efficiente: memorizza n + nnz indici o puntatori  • Adatto per y  ATx • risoluzione di problemi di programmazione lineare • Non adatto per y  Ax

7. Formati convenzionali: CSC

8. Il nuovo formato CSB • Consideriamo matrici sparse n×n con nnz elementi non nulli  • β = block-size parameter • valore ideale = circa √n • per semplicità si assume β potenza di 2 • CSB - Compressed Sparse Blocks • Partizionamento della matrice in blocchi quadrati di dimensione β × β • Numero di blocchi n2 / β 2 • Ordinamento Z-Morton interno ai blocchi  • Sostiene y  Ax e y  ATx

9. Il nuovo formato CSB

10. Il nuovo formato CSB

11. Il nuovo formato CSB

12. Il nuovo formato CSB

13. Il nuovo formato CSB

14. Prod. Matrice-Vettore con CSB CSB_SPMV (A, x, y) 1  for i0 to n/β−1 in parallel// riga di blocco 2     do Initialize a dynamic array Ri 3       Ri[0]0 // Array di indici per // gli ultimi blocchi dei chunk 4         count0 // Contatore nnz in un chunk 5         for j0 to n/β−2 6             do countcount+nnz(Ai,j) 7                 if count+nnz(Ai,j+1) > O(β) 8                     then// Fine chunk 9                     append j to Ri// Ultimo blocco del chunk 10                    count0 11        append n/β−1 to Ri 12        CSB_BLOCKROWV (A, i, Ri, x, y[(i∙β),…,((i+1)∙β)−1])

15. Prod. Matrice-Vettore con CSB Divisione in chunk -1

16. Prod. Matrice-Vettore con CSB Divisione in chunk -1

17. Prod. Matrice-Vettore con CSB Divisione in chunk -1

18. Prod. Matrice-Vettore con CSB Divisione in chunk -1

19. CSB_BLOCKROWV (A, i, R, x, y) 11  if R.length = 2 // Singolo chunk 12    then ℓR[0]+1 // Blocco più a sinistra nel chunk 13        rR[1] // Blocco più a destra nel chunk 14        if ℓ = r 15         then// Il chunk contiene un singolo blocco denso 16                startA.blk_ptr[ f(i,ℓ)] 17                endA.blk_ptr[ f(i,ℓ)+1]−1 18                CSB_BLOCKV (A, start, end, β, x, y) 19           else// Il chunk è sparso 20                multiply y(Ai,ℓ Ai,ℓ+1 … Ai,r)x serially 21        return // Se la riga di blocchi contiene più chunk 22  mid⌈R.length/2⌉−1 // Divide i chunk in due sottoinsiemi // Calcola il punto di split per il vettore x 23  xmid β·(R[mid]−R[0]) 24  Alloca un vettore z di cardinalità β, inizializzati a 0 25 in parallel 26      doCSB_BLOCKROWV(A, i, R[0…mid], x[0…xmid−1], y) 27      doCSB_BLOCKROWV(A, i, R[mid…R.length−1], x[xmid…x.length−1], z) 28  for k0 to β−1 29      do y[k]y[k]+z[k]

20. CSB_BLOCKROWV (A, i, R, x, y) 11  if R.length = 2 // Singolo chunk 12    then ℓR[0]+1 // Blocco più a sinistra nel chunk 13        rR[1] // Blocco più a destra nel chunk 14        if ℓ = r 15         then// Il chunk contiene un singolo blocco denso 16                startA.blk_ptr[ f(i,ℓ)] 17                endA.blk_ptr[ f(i,ℓ)+1]−1 18                CSB_BLOCKV (A, start, end, β, x, y) 19           else// Il chunk è sparso 20                multiply y(Ai,ℓAi,ℓ+1…Ai,r)x serially 21        return // Se la riga di blocchi contiene più chunk 22  mid⌈R.length/2⌉−1 // Divide i chunk in due sottoinsiemi // Calcola il punto di split per il vettore x 23  xmid β·(R[mid]−R[0]) 24  Alloca un vettore z di cardinalità β, inizializzati a 0 25 in parallel 26      doCSB_BLOCKROWV(A, i, R[0…mid], x[0…xmid−1], y) 27      doCSB_BLOCKROWV(A, i, R[mid…R.length−1],   x[xmid…x.length−1], z) 28  for k0 to β−1 29      do y[k]y[k]+z[k]

21. Prod. Matrice-Vettore con CSB y[β..2β-1] Split ricorsivo dei chunk

22. Prod. Matrice-Vettore con CSB y[β..2β-1] Split ricorsivo dei chunk

23. CSB_BLOCKV (A, start, end, dim, x, y) 28  if end−start ≤ O(dim) 29     then// Calcola la computazione seriale y←y+Mx 30        for kstart to end // A.val[start…end] è un blocco dim×dim 31           do y[A.row_ind[k]]  y[A.row_ind[k]] + A.val[k]·x[A.col_ind[k]] 32        return 33   // Ricorsione: divide il blocco M in 4 quadranti    34  binary search  start, start+1,…,end  per il più piccolo s2       tale che (A.row_ind[s2] & dim/2) ≠ 0 35  binary search start, start+1,…,s2−1  per il più piccolo s1       tale che (A.col_ind[s1] & dim/2) ≠ 0 36  binary search s2, s2+1,…,end  per il più piccolo s3       tale che (A.col_ind[s3] & dim/2) ≠ 0 37  in parallel 38    doCSB_BLOCKV (A, start, s1−1, dim/2, x, y)   // M00 39    doCSB_BLOCKV (A, s3, end, dim/2, x, y)       // M11 40  inparallel 41    doCSB_BLOCKV (A, s1, s2−1, dim/2, x, y)      // M01 42    doCSB_BLOCKV (A, s2, s3−1, dim/2, x, y)      // M10

24. Prod. Matrice-Vettore con CSB Decomposizione Z-Morton in 4 quadranti

25. Prod. Matrice-Vettore con CSB Decomposizione Z-Morton in 4 quadranti

26. Prod. Matrice-Vettore con CSB Decomposizione Z-Morton in 4 quadranti

27. Prod. Matrice-Vettore con CSB Decomposizione Z-Morton in 4 quadranti

28. Analisi di complessità • Al fine di valutare la complessità dell'algoritmo definiamo: • work: denotato con T1, rappresenta il tempo di esecuzione in una macchina monoprocessore monothread • span: denotato con T∞, rappresenta il tempo di esecuzione con un infinito numero di processi o thread • Viene definito grado di parallelismo il rapporto T1/T∞

29. Lemma 1 • Lemma 1: il formato CSR usa n∙log(nnz) + nnz∙log(n) bit di indici di supporto per una matrice n×n con nnz elementi non nulli. • Dimostrazione: per indicizzare x elementi sono necessari log(x) bit. Dal prodotto righe-colonne risultano n∙log(nnz) bit per il row_ptr e nnz∙log(n) bit per col_ind.

30. Lemma 1 (continua)

31. Lemma 2 • Lemma 2: Il formato CSB usa (n2/β2)∙log(nnz)+2∙nnz∙log(β) bit di indici di supporto per una matrice n×n con nnz elementi non nulli. • Dimostrazione: per ogni elemento in val, usiamo log(β) bit per rappresentare l'indice di riga e log(β) bit per rappresentare l'indice di colonna e richiede quindi nnz∙log(β) bit per ciascuno degli indici. Aggiungiamo lo spazio dato dall'array blk_ptr, ossia (n2/β2) ∙ log(nnz) bit.

32. Lemma 2 (continua)

33. Lemma 2 (continua)

34. Corollario 3 • Corollario 3: il formato CSB usa (n)∙log(nnz)+nnz∙log(n) bit di indici di supporto per una matrice n×n con nnz elementi non nulli, con β=√n.

35. Lemma 4 • Lemma 4: Su un blocco di dimensioni β×β, contenente r elementi non nulli, CSB_BlockV viene eseguito con work O(r) e span O(β). • Dimostrazione (span): • lo span relativo alla moltiplicazione di una matrice dim×dim può essere descritto da S(dim)=2∙S(dim/2)+O(log(dim))=O(dim). • 2∙S(dim/2): viene invocata 2 volte in parallelo la ricorsione su un singolo blocco di dim/2 • O(log(dim)): costo della ricerca binaria dei tre indici di split • caso base = O(dim):il caso base è seriale su O(dim) elementi ed è dominante sui casi ricorsivi

36. Lemma 4 (continua) • Dimostrazione (work): • Consideriamo l'albero di grado 4 generato dalle chiamate ricorsive della funzione CSB_BlockV; ogni nodo corrisponde alla computazione di un sottoblocco dim×dim, con dim=2h, e 0<h<log(β). • Se un nodo è una foglia, allora verifica il caso base ed ha sicuramente al più O(2h)=O(dim) elementi non nulli. Ne deriva quindi che il costo di computazione del nodo è O(r). • [..]

37. CSB_BLOCKV (A, start, end, dim, x, y) 28  if end−start ≤ O(dim) 29     then// Calcola la computazione seriale y←y+Mx 30        for kstart to end // A.val[start…end] è un blocco dim×dim 31           do y[A.row_ind[k]]  y[A.row_ind[k]] + A.val[k]·x[A.col_ind[k]] 32        return 33   // Ricorsione: divide il blocco M in 4 quadranti    34  binary search  start, start+1,…,end  per il più piccolo s2       tale che (A.row_ind[s2] & dim/2) ≠ 0 35  binary search start, start+1,…,s2−1  per il più piccolo s1       tale che (A.col_ind[s1] & dim/2) ≠ 0 36  binary search s2, s2+1,…,end  per il più piccolo s3       tale che (A.col_ind[s3] & dim/2) ≠ 0 37  in parallel 38    doCSB_BLOCKV (A, start, s1−1, dim/2, x, y)   // M00 39    doCSB_BLOCKV (A, s3, end, dim/2, x, y)       // M11 40  inparallel 41    doCSB_BLOCKV (A, s1, s2−1, dim/2, x, y)      // M01 42    doCSB_BLOCKV (A, s2, s3−1, dim/2, x, y)      // M10

38. Lemma 4 (continua)

39. Lemma 4 (continua)

40. Lemma 4 (continua)

41. Lemma 4 (continua) • [...] • Se un nodo è interno, allora ha almeno O(dim) elementi non nulli. • Il costo computazionale del nodo, senza considerare nodi figli, è pari a O(log(dim))=O(log(2h))=O(h), dovuto alla ricerca binaria. • I nodi generici di livello h sono al più O(r/dim), e quindi concorrono ad un lavoro complessivo per ogni livello pari a O(h∙r/dim). • sommando, per ogni h, nodi interni su tali livelli e nodi foglia, otteniamo O(r):

42. Lemma 5 • Questo lemma analizza la moltiplicazione fra una riga di blocchi e un vettore • Lemma 5: Su una riga di blocchi contenente n/β blocchi e r elementi non nulli, CSB_BlockrowV viene eseguito con work O(r) e span O(β∙log(n/β)).

43. Lemma 5 (continua) • Dimostrazione (work): • consideriamo la chiamata su una riga di blocco partizionata in C chunk e definiamo W(C) il lavoro (work) eseguito. La funzione inizializza un vettore z di O(β) elementi e richiama ricorsivamente se stessa due volte, sulla metà dell'input • Il lavoro è descritto dalla disequazioneW(C) ≤ 2∙W(⌈C/2⌉)+O(β)da cui deriva W(C)=O(C∙β+r), poiché si hanno C attivazioni di complessità O(β), più i casi base di complessità O(r), ossia la computazione seriale del prodotto. Il numero C di chunk è, al più, pari a O(r/β) nel caso in cui r=O(β)

44. Lemma 5 (continua) • Dimostrazione (span): • Lo span può essere descritto da S(C)=S(⌈C/2⌉)+O(β)=O(β∙log(C))+S(1) • Abbiamo che il caso base ha uno span pari O(β) sia nel caso della moltiplicazione seriale che in quella nella chiamata CSB_BlockV • Il caso base viene eseguito log(C) volte, con C ≤ n/β • Lo span complessivo è quindi O(β∙log(n/β))

45. Teorema 6 • Teorema 6: In una matrice n×n contenente r elementi non nulli, CSB_SpMV viene eseguito con un work O(r+n2/β2) e uno span di O(β∙lg(n/β))+n/β). • Dimostrazione: • CSB_SpMV ricostruisce i chunk e avvia la funzione CSB_BlockrowV. Il costo computazionale, per work e span deriva dal lemma precedente con l'aggiunta del costo necessario alla costruzione dei chunk. • [...]

46. Teorema 6 (continua) • Dimostrazione: • [...] • nel caso del work si aggiunge un costo pari a O(n2/β2) dovuto all'analisi di una singola riga di blocchi di costo O(n/β) per un costo totale di O(r+n2/β2) • nel caso dello span si aggiunge un costo pari a O(n/β) dato che è possibile parallelizzare l'operazione per ogni singola riga di blocchi, ottenendo un costo totale di O(β∙lg(n/β))+n/β)

47. Corollario 7 e 8 • Corollario 7: in una matrice n×n, contenente nnz≥ n valori non nulli, scegliendo β=O(√n), CSB_SpMV lavora con un work di O(nnz) e uno span di O(√n)∙log(n)) raggiungendo un parallelismo di almeno • O( nnz / (√n∙log(n)) ) • Corollario 8: in una matrice n×n, contenente nnz≥n valori non nulli, scegliendo β=O(√n), CSB_SpMV_T lavora con un work di O(nnz) e uno span di O(√n∙log(n)) raggiungendo un parallelismo di almeno • O( nnz / (√n∙log(n)) )

48. Lemma 9 • Lemma 9: In una matrice n×n, scegliendo β = O(√n), la serializzazione di CSB_SpMV richiede uno spazio di O(√n∙log(n)) non contando lo spazio occupato dalla matrice stessa • Dimostrazione:  • lo spazio complessivo utilizzato è dato da due overhead • il primo, R, è l'array dei chunk, che, per ogni riga, utilizza uno spazio O(n/β). Dato che β=√n, si ha che lo spazio complessivo utilizzato è O(√n). • il secondo, z, è il vettore temporaneo di dimensione pari a β. A causa della profondità della ricorsione, lo spazio utilizzato è O(β∙log(n))=O(√n∙log(n)) • la complessità finale è quindi O(√n∙log(n))+O(√n)=O(√n∙log(n))

49. Corollario 10 • Corollario 10: supponiamo un'esecuzione di CSB_SpMV per una matrice n×n con la scelta β=√n in un  work-stealing scheduler (preemptive round robin) con la proprietà busy-leaves. Allora l'esecuzione richiede uno spazio di O(P∙√n∙log(n)), con P pari al numero di processi utilizzati

50. Sperimentazione • Scelta valore di β • Performance media di Ax e ATx • Risultati reali sulla scalabilità degli algoritmi • matrici di medie dimensioni • matrici di grandi dimensioni