80
M´ ethodes directes pour la r´ esolution des syst` emes lin´ eaires
d’o` u on d´ eduit l’estimation voulue avec g(u) = nu/(1 − 2nu).
La strat´ egie du pivot, examin´ ee ` a la Section 3.5, permet de maˆ ıtriser la
taille des pivots et rend possible l’obtention d’estimations du type (3.38) pour
toute matrice.
3.3.3 Impl´ ementation de la factorisation LU
La matrice L ´ etant triangulaire inf´ erieure avec des 1 sur la diagonale et U ´ etant
triangulaire sup´ erieure, il est possible (et commode) de stocker directement
la factorisation LU dans l’emplacement m´ emoire occup´ e par la matrice A.
Plus pr´ ecis´ ement, U est stock´ ee dans la partie triangulaire sup´ erieure de A
(y compris la diagonale), et L occupe la partie triangulaire inf´ erieure stricte
(il est inutile de stocker les ´ el´ ements diagonaux de L puisqu’on sait a priori
qu’ils valent 1).
Le code MATLAB de l’algorithme est propos´ e dans le Programme 4. La
factorisation LU est stock´ ee directement ` a la place de la matrice A.
Programme 4 - lukji : Factorisation LU de la matrice A, version kji
function [A]=lukji(A)
% LUKJI Factorisation LU de la matrice A dans la version kji
% Y=LUKJI(A): U est stock´ e dans la partie triangulaire sup´ erieure
% de Y et L est stock´ e dans la partie triangulaire inf´ erieure
% stricte de Y.
[n,m]=size(A);
if n ˜= m, error(’Seulement les syst` emes carr´ es’); end
for k=1:n-1
if A(k,k)==0; error(’Pivot nul’); end
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n
i=[k+1:n]; A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
return
On appelle cette impl´ ementation de l’algorithme de factorisation version
kji, `
a cause de l’ordre dans lequel les boucles sont ex´ ecut´ ees. On l’appelle
´ egalement SAXP Y − kji car l’op´ eration de base de l’algorithme consiste `
a
effectuer le produit d’un scalaire par un vecteur puis une addition avec un
autre vecteur (SAXP Y est une formule consacr´ ee par l’usage ; elle provient
de “Scalaire A multipli´ e par vecteur X P lus Y ”).
La factorisation peut naturellement ˆ etre effectu´ ee dans un ordre diff´ erent.
Quand la boucle sur l’indice i pr´ ec` ede celle sur j, l’algorithme est dit orient´ e
ligne. Dans le cas contraire, on dit qu’il est orient´ e colonne. Comme d’habitude, cette terminologie provient du fait que la matrice est lue par lignes ou
par colonnes.
M´ ethodes directes pour la r´ esolution des syst` emes lin´ eaires
d’o` u on d´ eduit l’estimation voulue avec g(u) = nu/(1 − 2nu).
La strat´ egie du pivot, examin´ ee ` a la Section 3.5, permet de maˆ ıtriser la
taille des pivots et rend possible l’obtention d’estimations du type (3.38) pour
toute matrice.
3.3.3 Impl´ ementation de la factorisation LU
La matrice L ´ etant triangulaire inf´ erieure avec des 1 sur la diagonale et U ´ etant
triangulaire sup´ erieure, il est possible (et commode) de stocker directement
la factorisation LU dans l’emplacement m´ emoire occup´ e par la matrice A.
Plus pr´ ecis´ ement, U est stock´ ee dans la partie triangulaire sup´ erieure de A
(y compris la diagonale), et L occupe la partie triangulaire inf´ erieure stricte
(il est inutile de stocker les ´ el´ ements diagonaux de L puisqu’on sait a priori
qu’ils valent 1).
Le code MATLAB de l’algorithme est propos´ e dans le Programme 4. La
factorisation LU est stock´ ee directement ` a la place de la matrice A.
Programme 4 - lukji : Factorisation LU de la matrice A, version kji
function [A]=lukji(A)
% LUKJI Factorisation LU de la matrice A dans la version kji
% Y=LUKJI(A): U est stock´ e dans la partie triangulaire sup´ erieure
% de Y et L est stock´ e dans la partie triangulaire inf´ erieure
% stricte de Y.
[n,m]=size(A);
if n ˜= m, error(’Seulement les syst` emes carr´ es’); end
for k=1:n-1
if A(k,k)==0; error(’Pivot nul’); end
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
for j=k+1:n
i=[k+1:n]; A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
return
On appelle cette impl´ ementation de l’algorithme de factorisation version
kji, `
a cause de l’ordre dans lequel les boucles sont ex´ ecut´ ees. On l’appelle
´ egalement SAXP Y − kji car l’op´ eration de base de l’algorithme consiste `
a
effectuer le produit d’un scalaire par un vecteur puis une addition avec un
autre vecteur (SAXP Y est une formule consacr´ ee par l’usage ; elle provient
de “Scalaire A multipli´ e par vecteur X P lus Y ”).
La factorisation peut naturellement ˆ etre effectu´ ee dans un ordre diff´ erent.
Quand la boucle sur l’indice i pr´ ec` ede celle sur j, l’algorithme est dit orient´ e
ligne. Dans le cas contraire, on dit qu’il est orient´ e colonne. Comme d’habitude, cette terminologie provient du fait que la matrice est lue par lignes ou
par colonnes.