1.  Introduzione al calcolo parallelo in algebra lineare numerica


Definizione (un po' ironica, ma realistica) di super-calcolatore: " macchina con prestazioni pari ad 1/10 di quelle desiderate dagli utenti più esigenti " .
Dietro a questa definizione sta un fatto concreto: nonostante le impressionanti evoluzioni tecnologiche dell'hardware, le esigenze di calcolo risultano spesso superiori alla capacità espressa dal calcolatore più veloce disponibile.
Una risposta adeguata a questo problema viene data dal calcolo parallelo: più processori vengono messi a cooperare nella soluzione di un problema di calcolo. Facciamo qui riferimento alla configurazione cosiddetta " a memoria distribuita", in cui ogni processore ha una sua memoria a cui solo lui può accedere, ed un canale di comunicazione con gli altri processori (è la situazione che si ha, ad esempio, facendo cooperare più PC collegati mediante rete Ethernet);  l'altra configurazione di riferimento è quella cosiddetta " a memoria condivisa ", in cui più processori condividono l'accesso ad una medesima memoria (ed è la situazione che si ha, ad esempio, in una scheda bi-processore per PC).
La configurazione a memoria distribuita è più complessa da gestire e merita una trattazione più approfondita. Questa soluzione, richiede tre ingredienti principali:
Le reti di interconnessione esistenti sono molteplici e coprono un ampio spettro di prestazioni: si va dalla semplice rete Ethernet, con varie velocità, all'aggiunta di switch veloci del costo di qualche migliaio di Euro, fino alle reti ad alte prestazioni dei super-calcolatori del costo di qualche milione di Euro. In questa esercitazione si vedrà il peso della rete di interconnessione nelle prestazioni di calcolo e nella scelta dell'algoritmo più adatto. Se nella scorsa esercitazione abbiamo visto come sia necessario, data la grande velocità con la quale le CPU attuale effettuano operazioni di calcolo in virgola mobile, implementare gli algoritmi in modo che mantengano impegnata il più possibile la CPU in presenza di memoria gerarchica, qui (nella configurazione a memoria distribuita) il problema dello spostamento dei dati in memoria diventa ancora più importante.
Le librerie di comunicazione per il calcolo a memoria distribuita sono state sostanzialmente riunite in un unico standard, MPI  ( http://www-unix.mcs.anl.gov/mpi/ ).
Queste librerie permettono sostanzialmente di:
Gli algoritmi di calcolo parallelo di riferimento per l'algebra lineare numerica si trovano in www.netlib.org/scalapack  .

In generale, la buona riuscita di un algoritmo parallelo dipende da:




2.    Misura delle prestazioni nel calcolo parallelo

La misura di riferimento è generalmente il tempo di calcolo. Definiamo con:
T(p,n)    il tempo necessario a risolvere un problema di dimensione n  su  p  processori ;
T(1,n)    il tempo richiesto dal miglior algoritmo seriale ;
Speedup(p) = S(p) = T(1)/T(p)    dice quante volte si va più veloci con p processori, a pari n ;
Efficiency(p) = E(p) = S(p)/p    è l'efficienza, ed assume valori tra 0 e 1 .

In pratica, è raro riuscire ad avere S(p)=p per p elevato. E' più ragionevole richiedere che l'algoritmo sia scalabile, e cioè che l'efficienza sia limitata inferiormente al crescere di p (e quindi che aggiungendo processori si guadagni sempre qualcosa in velocità di esecuzione). Per capire se e fino a che punto un algoritmo è scalabile, è utile la legge di Amdahl : assumiamo, per semplicità, che il programma che implementa l'algoritmo contenga operazioni classificabili in due sole categorie: operazioni completamente parallelizzabili ed operazioni completamente seriali.  Detta f la frazione di programma completamente parallelizzabile  e detta s=1-f la frazione di programma completamente seriale, allora :
e dunque anche il numero di processori utilizzabile è limitato da 1/s   (parallelismo intrinseco).

Vediamo ora un modello (semplice, ma realistico) del processo di trasmissione dei dati da un processore all'altro, chiamato "LogP" . Si suppone che il tempo necessario a spedire un messaggio costituito da n pacchetti, sia il seguente:

2*o + L + (n-1)*g = (2*o + L - g) + n*g   =   alpha + n * beta

dove L è la latenza, cioè il tempo di percorrenza della rete, dal processore sorgente a quello destinazione;  o è il tempo che il processore deve spendere per spedire o ricevere un pacchetto di dati dalla rete, g è il minimo intervallo di tempo tra due spedizioni/ricezioni consecutive .
Queste grandezze vengono espresse relativamente all'unità di tempo che è pari al tempo necessario ad un processore ad eseguire un'operazione in virgola mobile.
Nella pratica, si ha che  alpha » 1 e beta » 1 e quindi il tempo necessario a trasferire dati tra i processori è molto più grande (anche centinaia di migliaia di volte) del tempo necessario ad eseguire un'operazione in virgola mobile. Inoltre,  alpha » beta  e quindi è meglio fare pochi messaggi lunghi piuttosto che tanti e brevi . In generale, i tempi della comunicazione sono dovuti in gran parte al software, non all'hardware.




3.   Esempio: prodotto di matrici usando differenti reti di connessione

Ora che abbiamo a disposizione un modello per il calcolo dei tempi di trasmissione, lo utilizziamo sull'esempio del prodotto di matrici, C = A*B, utilizzando differenti reti di interconnessione e differenti distribuzioni dei dati.

a.    bus senza broadcast (es. Ethernet) con suddivisione a blocchi di colonne:  ogni processore possiede sempre un blocco di colonne [n x n/p] di ognuna delle matrici A, B e C; dato che un blocco di colonne di C risulta dalla moltiplicazione dell'intera matrice A per il corrispondente blocco di colonne di B (vedi Algoritmo BLAS2 della Esercitazione n.2 del corso di C.N.), è necessario che ogni processore riceva dagli altri, in sequenza, i p-1 blocchi di colonne di A diversi da quello posseduto all'inizio. Il costo di spedizione sulla rete di un blocco di colonne è pari a   t_c = alpha + n * (n/p) * beta   .   Relativamente a ciascun blocco di A, il processore esegue operazioni per  t_a = 2 * n^3 / p^2   .
Dato che, per p fissato,  t_c = O(n^2)  e  t_a = O(n^3)  ,  si ha che da un certo  n  in poi, in genere elevato ,  t_a domina su t_c e in tal caso i processori non rimangono fermi in attesa dei dati.

Algoritmo:    C_Pmio = C_Pmio + A_Pmio * B_Pmio_Pmio ;
                    for i=0:Pmio-1        aspetta di ricevere A_i  e poi calcola  C_Pmio = C_Pmio + A_i * B_i_Pmio ;
                    for i=0:p-1 eccetto Pmio        spedisce A_Pmio  al processore "i"
                    for i=Pmio+1:p-1        aspetta di ricevere A_i  e poi calcola  C_Pmio = C_Pmio + A_i * B_i_Pmio ;

che produce il seguente diagramma temporale:
bus_senza_broadcast

( versione testuale del diagramma in figura :    ESE_CN_Informatica_3_diagramma_a.txt   )

da cui si ha:        T(p) = p * (p-1) * t_c  +  3 * t_a    e dunque, essendo  T(1) = t_a * p^2  ,  E(p) = 1 / (3/p + (p-1)*t_c/t_a)
e dunque l'algoritmo è scalabile (cioè l'efficienza non tende a zero per p crescente) solo se  t_c ≈ t_a / (p-2)  ,  che equivale alla condizione di pieno utilizzo della rete (se t_c è inferiore, si hanno buchi di inattività della rete, se è superiore, si hanno maggiori tempi di inattività dei processori). Essendo in generale la rete una risorsa relativamente lenta, cercare di saturare la rete è una buona regola generale di programmazione numerica parallela.





b.    bus con broadcast e con suddivisione a blocchi di colonne:  analogo al caso precedente.
Algoritmo:      C_Pmio = C_Pmio + A_Pmio * B_Pmio_Pmio ;
                        for  i=0:p-1  
                            if   i==Pmio    broadcast  A_Pmio
                            else  aspetta di ricevere A_i  e poi calcola  C_Pmio = C_Pmio + A_i * B_i_Pmio ;

che produce il seguente diagramma temporale:
bus_con_broadcast

( versione testuale del diagramma in figura :    ESE_CN_Informatica_3_diagramma_b.txt   )

da cui si ha:        T(p) = p * t_c  +  (p+1) * t_a    e dunque, essendo  T(1) = t_a * p^2  ,  E(p) = 1 / ((p+1)/p + t_c/t_a)
e dunque, facendo dominare t_a  su  t_c ,  si può ottenere un'efficienza sensibilmente più grande per valori di  n  sensibilmente inferiori al caso precedente.





c.    connessione ad anello (ring) e con suddivisione a blocchi di colonne:  la connessione ad anello la si ottiene collegando ogni processore ad un suo confinante sinistro (Pmio-1) e ad un confinante destro (Pmio+1) in modo da costruire un percorso chiuso; fisicamente la si può realizzare con un insieme di PC ognuno con due schede di rete .

Algoritmo:    copia A_Pmio su T ;
                    C_Pmio = C_Pmio + A_Pmio * B_Pmio_Pmio ;
                    for  i = 1  to  p-1
                           spedisce T  al processore (Pmio+1) mod p ;
                           aspetta di ricevere T  dal processore "rem(Pmio-1, p)"
                           C_Pmio = C_Pmio + T * B_rem(Pmio-1,p)_Pmio ;
                           
da cui si ha :    T(p) = (p-1) * t_c  +  p * t_a  e dunque    E(p) = 1 / ( 1 + (p-1)/p * t_c / t_a)
che è un risultato analogo al precedente, ma con una connessione di costo nettamente inferiore.


A questo punto, notiamo che il risultato raggiunto è ottimale perchè, per come sono distribuiti i dati, è sempre e comunque necessario che (p-1) blocchi di A arrivino ad ogni processore. Quindi, per migliorare è necessario poter distribuire e trasmettere i dati in maniera differente, come nella modalità seguente:


d.    connessione a griglia bi-dimensionale, con suddivisione a blocchi quadrati:  si suppone che le matrici A, B e C siano partizionate in blocchi quadrati di dimensione [n/s x n/s],  con  p = s^2 , e che ogni processore possegga un blocco di A, uno di B ed uno di C . I processori sono disposti ai nodi di una griglia quadrata, di lato lungo s ed in cui ogni riga ed ogni colonna sono chiuse ad anello.
L'algoritmo utilizzato è dovuto a Cannon, e si basa sul seguente riordino delle operazioni:     invece di eseguire  
                                                                                                                  s-1
C_i_j   =   Σ    ( A_i_k * B_k_j )
                                                                                                               k = 0
si  esegue  :
                                                                                                 s-1
C_i_j   =   Σ    ( A_i_rem(i+j+k,s) * B_rem(i+j+k,s)_j )
                                                                                              k = 0

in pratica, viene cambiato l'ordine con cui viene eseguita la somma. Questo può essere realizzato tramite uno spostamento iniziale delle matrici, da (es. per p=9):

                            A_0_0        A_0_1        A_0_2                                    A_0_0        A_0_1        A_0_2       
                            A_1_0        A_1_1        A_1_2                a :                A_1_1        A_1_2        A_1_0       
                            A_2_0        A_2_1        A_2_2                                    A_2_2        A_2_0        A_2_1       

                            B_0_0        B_0_1        B_0_2                                    B_0_0        B_1_1        B_2_2
                            B_1_0        B_1_1        B_1_2                                    B_1_0        B_2_1        B_0_2
                            B_2_0        B_2_1        B_2_2                                    B_2_0        B_0_1        B_1_2

e poi da un ciclo iterativo (s iterazioni) composto da due operazioni (che possono essere svolte anche in sovrapposizione):
1)  ogni processore calcola il prodotto dei blocchi di A e di B in suo possesso ed aggiorna il suo blocco di C ;
2)  ogni processore spedisce a sinistra il suo blocco di A e verso l'alto il suo blocco di B, e contemporaneamente riceve da destra e dal basso i nuovi blocchi (notare che questa operazione sature la rete di interconnessione di ogni processore, formata da quattro segmenti: alto, basso , sinistro e destro; in teoria la si potrebbe realizzare usando PC con quattro schede di rete).
Si ottiene :    T(p) = sqrt(p) * t_c  +  p * t_a  e dunque    E(p) = 1 / ( 1 + 1/sqrt(p) * t_c / t_a)
che è sensibilmente migliore dei precedenti .