Applicazione della decomposta incompleta

Mediante l'algoritmo sviluppato nel precedente paragrafo non si calcola esplicitamente la matrice di precondizionamento $K^{-1}$ ma solamente il fattore incompleto $\tilde{U}$. È quindi necessario sviluppare un algoritmo ad hoc che permetta di calcolare il generico prodotto $K^{-1} \bv$ senza generare esplicitamente $K^{-1}$.

Sia $\bw$ il vettore in $\Re^n$ risultato del prodotto $K^{-1} \bv$. Per definizione di $K^{-1}$ si ha:

$$\bw = \left( \tilde{U}^T \tilde {U} \right)^{-1} \bv$$ (27)

cioè, premoltiplicando per $K$ ambo i membri:

$$\left( \tilde{U}^T \tilde{U} \right) \bw = \bv$$ (28)

Il calcolo di $\bw$ viene quindi ricondotto alla soluzione di un sistema lineare la cui matrice è $K$. Poiché $K$ è fattorizzabile nel prodotto di due matrici triangolari, il sistema (28) è risolvibile tramite sostituzioni in avanti e all'indietro. Posto:

$$\tilde{U} \bw = \bz$$ (29)

la (28) diventa:

$$\tilde{U}^T \bz = \bv$$ (30)

che si può facilmente risolvere con sostituzioni in avanti (Figura 1). Iniziando dalla prima componente, ricavata in modo immediato come:

$$z_1 = \frac{v_1}{u_{11}}$$ (31)

si ottiene con semplici calcoli la formula ricorrente:

$$z_i = \frac{1}{u_{ii}} \left( v_i - \displaystyle{\sum_{j=1}^{i-1}} u_{ji} z_j \right) \;\;\;\;\;\;\; i=2,\ldots,n$$ (32)

Figura 1: Schema della soluzione del sistema con sostituzioni in avanti.

Come osservato per il prodotto matrice-vettore e per la determinazione del fattore incompleto di Cholesky, la sommatoria contenuta in (32) non è banale da implementare in modo efficiente memorizzando le matrici secondo il sistema CRS. A tal proposito, risulta conveniente definire un vettore di accumulazione $\bs$ nelle cui componenti viene aggiornato il prodotto contenuto nella sommatoria dell'equazione (32) procedendo per righe di $\tilde{U}$. L'algoritmo per il calcolo del vettore $\bz$ può pertanto essere scritto come:

001 Per $j=1,n$
002 azzero $s_j$
003 Fine Per
004 $z_1:=v_1/$PREC$(1)$
005 Per $i=2,n$
006 $k:=$ IA$(i)$
006 Per $k1=$ IA$(i-1)+1$, IA$(i)-1$
007 $j:=$ JA$(k1)$
008 $s_j:=s_j \; +$ PREC $(k1) \cdot z_{i-1}$
009 Fine Per
010 $z_i:=\left(v_i-s_i\right)/$PREC$(k)$
011 Fine Per

Ottenuto il vettore $\bz$ si può infine calcolare $\bw$ risolvendo il sistema (29) tramite sostituzioni all'indietro (Figura 2). La formula ricorrente si ricava in modo del tutto analogo a quanto fatto nelle equazioni (31) e (32) partendo in questo caso dalla componente $n$-esima:

$$w_n = \frac{z_n}{u_{nn}}$$ (33)

$$w_i = \frac{1}{u_{ii}} \left( z_i - \displaystyle{\sum_{j=i+1}^{n}} u_{ij} w_j \right) \;\;\;\;\;\;\; i=n-1,\ldots,1$$ (34)

Figura 2: Schema della soluzione del sistema con sostituzioni all'indietro.

L'implementazione della (2) risulta stavolta molto semplice anche memorizzando la matrice $\tilde{U}$ in forma compatta. Sempre utilizzando il vettore di accumulo $\bs$, l'algoritmo corrispondente può essere scritto nel modo seguente:

001 $w_n:=z_n/$PREC$(nt)$
002 Per $i=n-1,1$ con passo -1
003 azzero $s_i$
004 $k:=$ IA$(i)$
005 Per $k1=$ IA$(i)+1$, IA$(i+1)-1$
006 $j:=$ JA$(k1)$
007 $s_i:=s_i \; +$ PREC $(k1) \cdot w_j$
008 Fine Per
009 $w_i:=\left(z_i-s_i\right)/$PREC$(k)$
010 Fine Per

Per un confronto, viene messa a disposizione dello studente la subroutine LSOLVE che implementa gli algoritmi soprariportati.