Decomposizione LU

1
Decomposizione LU
Consideriamo un sistema di n equazioni con n incognite:
a0,0 x0 + a0,1 x1 + · · · + a0,n−1 xn−1
a1,0 x0 + a1,1 x2 + · · · + a1,n−1 xn−1
=
=
..
.
b0
b1
an−1,0 x0 + an−1,1 x1 + · · · + an−1,n−1 xn−1
=
bn−1
Ponendo

a0,0
 ..
A= .
an−1,0
···

a0,n−1

..
,
.
· · · an−1,n−1


x0


x =  ...  ,
xn−1
Il sistema pu`
o essere scritto in forma matriciale:
Ax = b.


b0


b =  ... 
bn−1
(1)
La decomposizione LU `e una metodo che, a partire da una matrice A ∈ Rn×n genera due matrici
L, U ∈ Rn×n , con L triangolare inferiore a elementi diagonali unitari e U triangolare superiore,
tali che A = LU .
In questo modo `e possibile risolvere il sistema lineare Ax = b risolvendo due sistemi lineari
triangolari:
(
Ly = b
Ux = y
Determinate le matrici L e U , prima si risolve il sistema con matrice dei coefficienti L, termine
noto b e incognita y attraverso l’algoritmo di eliminazione in avanti; poi si risolve il sistema che
ha come matrice U , termine noto il vettore y appena determinato e incognita il vettore x con
l’algoritmo di sostitituzione all’indietro.
L’algoritmo della decomposizione LU si basa sull’eliminazione di Gauss: il principio `e quello
di eliminare progressivamente delle variabili dal sistema, finch`e non rimane una sola variabile.
Presa ad esempio la i-esima equazione nel sistema, possiamo eliminare un’incognita da un’altra
equazione j, moltiplicando l’equazione i per il termine costante −aji /aii e aggiungendo l’equazione risultante all’equazione j. L’obiettivo `e quello di azzerare tutti i termini che stanno
nel triangolo inferiore della matrice, procedendo per colonne, ed ottenendo cos`ı la matrice triangolare superiore U . I coefficienti aji /aii , per cui vengono moltiplicate le equazioni, andranno
a definire il triangolo inferiore della matrice L.
La versione sequenziale dell’algoritmo si pu`
o scrivere in questo modo:
FOR ii = 0,..,n-1
FOR jj = ii+1,...,n-1
m = a(jj,ii) / a(ii,ii)
FOR kk = ii+1, ...,n-1
a(jj,kk) -= m * a(ii,kk)
END FOR
a(jj,ii) = m
END FOR
END FOR
2
n3
. L’accesso ai dati nel ciclo pi`
u interno
3
avviene per righe, di conseguenza si ha un buono sfruttamento della cache.
Nel caso considerato, l’algoritmo `e descritto in forma diagonale, utilizzabile solo per particolari tipi di matrici. Un algoritmo pi`
u generale invece attua la strategia del pivoting. Nell’esercitazione prenderemo in considerazione solo matric strettamente diagonali dominanti, per le
quali cio`e vale:
Il costo computazionale dell’algoritmo `e O
|ai,i | >
n−1
X
|ai,j |
j=0,j6=i
In questo caso, `e sufficiente l’algoritmo di Gauss con strategia diagonale.
Vediamo ora come `e possibile parallelizzare questo algoritmo. Supponiamo di avere un
processore per ogni riga della matrice, ovvero p = n. Ogni processore memorizza gli elementi
della riga a lui associata, ovvero il processore Pi ha gli elementi della riga ai .
L’algoritmo parallelo avr`
a n − 1 iterazioni: nella prima iterazione, il processore P0 invia
a tutti gli altri, tramite una comunicazione broadcast, gli elementi della propria riga. Tutti
gli altri processori Pi : i > 0 potranno cos`ı azzerare il primo elemento della propria riga. Nella
seconda iterazione, il processore P1 invia a tutti gli altri, tramite una comunicazione broadcast,
gli elementi della propria riga, eccetto il primo elemento perch`e `e gi`a stato azzerato. Tutti
gli altri processori Pi : i > 1 potranno cos`ı azzerare il secondo elemento della propria riga. In
generale, alla j-esima iterazione, il processore Pj comunica le sue rimanenti n − j componenti
ai processori Pi : i > j, i quali azzereranno la loro (j + 1)-esima componente.
Si vede chiaramente come molti processori sono idle, ovvero non compiono nessuna operazione. Ad esempio, il processore P0 , dopo aver comunicato la propria riga agli altri processori,
non ha pi`
u compiti da eseguire per tutto il resto dell’algoritmo.
Il carico di lavoro `e quindi mal bilanciato e lo speedup di questo algoritmo non potr`
a essere
buono. Ci si poteva aspettare un risultato simile, in quanto siamo partiti dall’assunto di avere
tanti processori quante sono le righe della matrice. Il numero di processori sarebbe troppo
elevato per ottenere buoni speedup.
Analizziamo ora un caso pi`
u realistico, ovvero quello in cui si hanno molti meno processori
rispetto alle righe della matrice. Questo `e un caso pi`
u reale, in quanto generalmente le matrici da fattorizzare hanno dimensioni nell’ordine delle migliaia, mentre il numero di processori
disponibili `e nell’ordine della decina.
Se distribuissimo la matrice A a blocchi per righe non risolveremmo il problema del cattivo
bilanciamento del carico. Il processore P0 avrebbe pochissime operazioni da svolgere, dato che
dovrebbe azzerare gli elementi del triangolo superiore. L’ultimo processore dovrebbe al contrario
operare su quasi tutti i suoi elementi.
Per bilanciare il carico di lavoro, `e allora utile utilizzare una distribuzione ciclica delle righe,
ovvero il processore Pj possiede le righe {ai : i = j mod p}, dove p `e il numero totale di
processori (vedi Figura ).
Con questa distribuzione, l’algoritmo di decomposizione LU parallelo si pu`
o schematizzare
nel seguente modo:
1. Il processore P0 comunica agli altri processori le righe della matrice da fattorizzare, secondo
la distribuzione ciclica. La riga i sar`a inviata al processore i mod p.
2. Consideriamo la riga i: il processore i mod p, ovvero il processore che detiene la riga
i-esima, comunica con una broadcast la riga a tutti i processori.
3
Figura 1:
3. Sia Rq = {q + kp : k ∈ N, q + kp < n} l’insieme delle righe possedute dal processore
q-esimo. Per ciascuna riga in Rq , con indice maggiore di i, annulla l’i-esimo elemento,
sfruttando la riga ricevuta al passo 2.
4. Memorizzare, sempre nella matrice A, i coefficienti usati per annullare gli elementi.
5. Ripetere dal passo 2, per i = 0, . . . , n − 2.
6. Ogni processore comunica a P0 le proprie righe della matrice A (che ora conterr`a gli
elementi di L e U .
` possibile osservare che in questo algoritmo `e necessario operare un grande numero di coE
municazioni: occorrono n comunicazioni per il passo 1, altre n per il passo 2 (una comunicazione
ogni iterazione) e n per il passo 6. Quindi, in totale, il numero di comunicazioni aumenta all’aumentare della dimensione del problema. Questa `e una caratteristica atipica poich´e solitamente
il numero di comunicazione `e proporzionale al numero di processori e non alla dimensione del
problema (aumenta la quantit`
a di dati trasmessi, ma non il numero di comunicazioni). Inoltre, la
distribuzione ciclica produce un bilanciamento perfetto del carico di lavoro solo asintoticamente,
quindi pi`
u sono piccole le dimensioni della matrice A, pi`
u si sar`a lontani da un bilanciamento
ottimale. Considerati questi punti, non sar`a possibile osservare speedup ottimali per questo
algoritmo.

Download Report