132
4 M´ ethodes it´ eratives pour la r´ esolution des syst` emes lin´ eaires
2. Factorisation LU incompl` ete (ILU en abr´ eg´ e) et Factorisation de Cholesky
incompl` ete (IC en abr´ eg´ e).
Une factorisation LU incompl` ete de A consiste ` a calculer une matrice
triangulaire inf´ erieure L in et une matrice triangulaire sup´ erieure U in (approximations des matrices exactes L et U de la factorisation LU de A),
telles que le r´ esidu R = A−L in U in poss` ede des propri´ et´ es donn´ ees, comme
celle d’avoir certains coefficients nuls. L’objectif est d’utiliser les matrices
L in , U in comme pr´ econditionneur dans (4.24), en posant P = L in U in .
Nous supposons dans la suite que la factorisation de la matrice A peut
ˆ etre effectu´ ee sans changement de pivot.
L’id´ ee de base de la factorisation incompl` ete consiste ` a imposer ` a la matrice approch´ ee L in (resp. U in ) d’avoir la mˆ eme structure creuse que
la partie inf´ erieure (resp. sup´ erieure) de A. Un algorithme g´ en´ eral pour
construire une factorisation incompl` ete est d’effectuer une ´ elimination de
Gauss comme suit : `
a chaque ´ etape k, calculer m ik = a
(k)
ik /a
(k)
kk seulement si a ik = 0 pour i = k + 1, . . ., n. Calculer alors a
(k+1)
ij
pour
j = k + 1, . . . , n seulement si a ij = 0. Cet algorithme est impl´ ement´ e
dans le Programme 17 o` u la matrice L in (resp. U in ) est progressivement
´ ecrite ` a la place de la partie inf´ erieure de A (resp. sup´ erieure).
Programme 17 - basicILU : Factorisation LU incompl` ete (ILU)
function [A] = basicILU(A)
%BASICILU Factorisation LU incompl` ete.
% Y=BASICILU(A): U est stock´ ee dans la partie triangulaire sup´ erieure
% de Y et L dans la partie triangulaire inf´ erieure stricte de Y.
% Les matrices L et U ont la mˆ eme structure creuse que la matrice A
[n,m]=size(A);
if n ˜= m, error(’Seulement pour les matrices carr´ ees’); end
for k=1:n-1
for i=k+1:n,
if A(i,k) ˜= 0
if A(k,k) == 0, error(’Pivot nul’); end
A(i,k)=A(i,k)/A(k,k);
for j=k+1:n
if A(i,j) ˜= 0
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
end
end
return
Remarquer que le fait d’avoir la mˆ eme structure creuse pour L in (resp.
U in ) que pour la partie inf´ erieure (resp. sup´ erieure) de A, n’implique pas
4 M´ ethodes it´ eratives pour la r´ esolution des syst` emes lin´ eaires
2. Factorisation LU incompl` ete (ILU en abr´ eg´ e) et Factorisation de Cholesky
incompl` ete (IC en abr´ eg´ e).
Une factorisation LU incompl` ete de A consiste ` a calculer une matrice
triangulaire inf´ erieure L in et une matrice triangulaire sup´ erieure U in (approximations des matrices exactes L et U de la factorisation LU de A),
telles que le r´ esidu R = A−L in U in poss` ede des propri´ et´ es donn´ ees, comme
celle d’avoir certains coefficients nuls. L’objectif est d’utiliser les matrices
L in , U in comme pr´ econditionneur dans (4.24), en posant P = L in U in .
Nous supposons dans la suite que la factorisation de la matrice A peut
ˆ etre effectu´ ee sans changement de pivot.
L’id´ ee de base de la factorisation incompl` ete consiste ` a imposer ` a la matrice approch´ ee L in (resp. U in ) d’avoir la mˆ eme structure creuse que
la partie inf´ erieure (resp. sup´ erieure) de A. Un algorithme g´ en´ eral pour
construire une factorisation incompl` ete est d’effectuer une ´ elimination de
Gauss comme suit : `
a chaque ´ etape k, calculer m ik = a
(k)
ik /a
(k)
kk seulement si a ik = 0 pour i = k + 1, . . ., n. Calculer alors a
(k+1)
ij
pour
j = k + 1, . . . , n seulement si a ij = 0. Cet algorithme est impl´ ement´ e
dans le Programme 17 o` u la matrice L in (resp. U in ) est progressivement
´ ecrite ` a la place de la partie inf´ erieure de A (resp. sup´ erieure).
Programme 17 - basicILU : Factorisation LU incompl` ete (ILU)
function [A] = basicILU(A)
%BASICILU Factorisation LU incompl` ete.
% Y=BASICILU(A): U est stock´ ee dans la partie triangulaire sup´ erieure
% de Y et L dans la partie triangulaire inf´ erieure stricte de Y.
% Les matrices L et U ont la mˆ eme structure creuse que la matrice A
[n,m]=size(A);
if n ˜= m, error(’Seulement pour les matrices carr´ ees’); end
for k=1:n-1
for i=k+1:n,
if A(i,k) ˜= 0
if A(k,k) == 0, error(’Pivot nul’); end
A(i,k)=A(i,k)/A(k,k);
for j=k+1:n
if A(i,j) ˜= 0
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
end
end
return
Remarquer que le fait d’avoir la mˆ eme structure creuse pour L in (resp.
U in ) que pour la partie inf´ erieure (resp. sup´ erieure) de A, n’implique pas