Fattore incompleto di Cholesky secondo Kershaw

Il calcolo del fattore incompleto di Cholesky si basa sulla fattorizzazione triangolare di matrici simmetriche che viene di seguito brevemente richiamata. Poichè con la memorizzazione compatta di $A$ vengono conservati i soli elementi appartenenti alla triangolare alta, conviene riscrivere la (21) come:

$$K^{-1} = \left( \tilde{U}^T \tilde{U} \right)^{-1}$$ (22)

dove $\tilde{U}=\tilde{L}^T$.

Supponiamo che $A$ sia una matrice 3$\times$3 piena di cui svolgiamo per esteso la fattorizzazione:

$$\left[ \begin{array}{ccc} u_{11} & 0 & 0 u_{12} & u_{22} & 0 \... ...13} a_{12} & a_{22} & a_{23} a_{13} & a_{23} & a_{33} \end{array} \right]$$ (23)

Si noti che nella (23) si è sfruttata la simmetria di $A$ ed il fatto che i soli elementi $a_{ij}$ memorizzati sono quelli per cui $j \geq i$.

Sviluppiamo il prodotto (23) procedendo secondo la successione delle righe di $A$. Per la prima riga vale:

$$a_{11} = u_{11}^2 \Rightarrow u_{11} = \sqrt{a_{11}}$$

$$a_{12} = u_{11} u_{12} \Rightarrow u_{12} = \frac{a_{12}}{u_{11}}$$ (24)

$$a_{13} = u_{11} u_{13} \Rightarrow u_{13} = \frac{a_{13}}{u_{11}}$$

Si può osservare che i coefficienti sottodiagonali di $A$ vengono di volta in volta già utilizzati nel calcolo delle righe superiori a quella corrente. Si procede, pertanto, con le righe 2 e 3 considerando i soli termini relativi alla triangolare alta di $A$:

$$a_{22} = u_{12}^2 + u_{22}^2 \Rightarrow u_{22} = \sqrt{a_{22} - u_{12}^2}$$

$$a_{23} = u_{12} u_{13} + u_{22} u_{23} \Rightarrow u_{23} = \frac{1} {u_{22}} \left( a_{23} - u_{12} u_{13} \right)$$ (25)

$$a_{33} = u_{13}^2 + u_{23}^2 + u_{33}^2 \Rightarrow u_{33} = \sqrt{a_{33} - u_{13}^2 - u_{23}^2}$$

Dalle (24) e (25) si può facilmente generalizzare l'algoritmo di calcolo del fattore completo di Cholesky per una matrice di ordine $n$:

001 $u_{11}:=\sqrt{a_{11}}$
002 Per $j=2,n$
003 $u_{1j}:=a_{1j}/u_{11}$
004 Fine Per
005 Per $i=2,n$
006 $u_{ii}:=\sqrt{a_{ii}-\sum_{l=1}^{i-1} u_{li}^2}$
007 Per $j=i+1,n$
008 $u_{ij}:=\left(a_{ij}-\sum_{l=1}^{i-1} u_{li}\;u_{lj}\right)/ u_{ii}$
009 Fine Per
010 Fine Per

Il passaggio concettuale dal fattore completo di Cholesky a quello incompleto secondo Kershaw risulta a questo punto banale, in quanto all'algoritmo precedente è sufficiente aggiungere l'assegnazione:

$$a_{ij}=0 \Rightarrow u_{ij}:=0$$ (26)

Poiché $A$ è memorizzata in modo compatto, e così anche $\tilde{U}$, si deduce immediatamente che, in virtù dell'assegnazione (26), i vettori IA e JA descrivono anche la topologia di $\tilde{U}$ per la quale sarà sufficiente adottare un vettore PREC contenente i coefficienti non nulli memorizzati sequenzialmente per righe.

Il nuovo algoritmo di calcolo del fattore incompleto di Cholesky secondo Kershaw con la memorizzazione in formato CRS risulta pertanto:

001 PREC $(1):=\sqrt{\mbox{\tt SYSMAT}(1)}$
002 Per $k=$ IA(1)+1, IA(2)-1
003 PREC$(k):=$ SYSMAT$(k)/$PREC(1)
004 Fine Per
005 Per $i=2,n$
006 $k:=$ IA$(i)$
007 $l:=$ ogni $k<{\mbox{\tt IA}(i)}$ per cui JA$(k)=i$
008 PREC $(k):=\sqrt{\mbox{\tt SYSMAT}(k)-\sum_l \mbox{\tt PREC}(l)^2}$
009 Per $k1=$ IA$(i)+1,$ IA$(i+1)-1$
010 $j:=$ JA$(k1)$
011 $l1:=$ ogni $k<{\mbox{\tt IA}(i)}$ per cui JA$(k)=i$
012 $l2:=$ ogni $k<{\mbox{\tt IA}(i)}$ per cui JA$(k)=j$
013 PREC $(k1):=\left(\mbox{\tt SYSMAT}(k1)-\sum_{l1,l2} \mbox{\tt PREC}(l1) \cdot \mbox{\tt PREC}(l2) \right)/\mbox{\tt PREC}(k)$
014 Fine Per
015 Fine Per

Va osservato, tuttavia, che l'implementazione efficiente delle sommatorie in riga 8 e 13 sugli indici $l$, $l1$ ed $l2$ è in realtà non banale e va preferibilmente effettuata calcolando separatamente il contributo del coefficiente relativo al $k$ corrente in modo da evitare le dispendiose ricerche previste in riga 7, 11 e 12. Tale implementazione è effettuata nella subroutine KERSH a disposizione dello studente.

Si deve infine ricordare che l'uso generalizzato delle (24) e (25) comporta l'estrazione di radici quadrate il cui argomento, nel caso di una fattorizzazione incompleta, non è più garantito essere positivo. In questo caso l'elemento diagonale in questione verrà posto pari ad un numero positivo arbitrario (ad esempio, l'ultimo coefficiente diagonale di $\tilde{U}$ non nullo). Si dimostra comunque teoricamente che se la matrice $A$ è di tipo M, cioè diagonalmente dominante con elementi extra-diagonali di segno opposto rispetto a quelli diagonali, come ad esempio si verifica nella discretizzazione agli elementi finiti dell'equazione ellittica di Laplace, tutte le radici da calcolare nel fattore incompleto esistono nel campo reale.