Initialisation

Nous aurons besoin de n pour la taille de la matrice, et de perm pour la permutation des indices des colonnes. Nous aurons aussi besoins de quelques autres variables de type integer et real que nous déclarons maintenant. La construction (/ … /) permet d’écrire un vecteur explicite, avec des coordonnées séparées par des virgules. La constrution (…,i=i_0,i_n) crée une suite de nombres séparés par une virgule qui se calculent avec ce qui remplacent les trois petits points à l’aide d’un indice comme dans une boucle do…end do : ici l’indice est i et il varie de i_0 à i_n. Le reste de l’algorithme est une boucle répétée n 1 fois (formellement, on verra qu’il y a n étapes, mais la dernière ne contient aucune instruction).
hDécompositon A = LUP 1i+
(9) /1
integer, dimension(size(m,2)) :: perm
integer :: n,j,k,ind_piv
real :: r,piv
n=size(m,1)
! m doit etre carree
perm=(/(k,k=1,n)/) ! initialise la permutation a l’identite : perm(k)=k do k=1,n-1
hpivotage dans la ligne 3i
hcalculs 4i
end do
end subroutine aelup

Pivotage partiel

En vue minimiser les erreurs de calculs liées à la représentation en machine des nombres réels, on
choisit à l’étape k de chercher dans la ligne k quelle colonne porte le plus grand élément (en valeur
absolue). Ensuite, plutôt que d’échanger les contenus des blocs mémoire de la matrice, on échange
les éléments de perm. Ainsi perm(k) est l’indice de la colonne qui porterait le numéro k si l’on avait
effectué les permutations de colonnes. Le coefficient M(i; j) d’indices (i; j) de la matrice M permutée
est donc m(i,perm(j)).
hpivotage dans la ligne 3i
(2)
ind_piv=k
piv=abs(m(k,perm(ind_piv)))
do j=k+1,n
r=abs(m(k,perm(j)))
if (r>piv) then
ind_piv=j
piv=r
end if
end do
trans(k)=ind_piv ! enregistre l’indice du j-eme pivot
if (ind_piv/=k) then
! echange perm(k) et perm(ind_piv)
j=perm(k);perm(k)=perm(ind_piv);perm(ind_piv)=j end if

Élimination de Gauss

Pour justifier l’algorithme (sans pivotage pour simplifier), donnons un invariant. On va fabriquer une suite de trois matrices carrées (Ak; Lk; Uk)nk=1+1 avec pour tout k l’égalité Ak + LkUk = A ; au début A1 = A, L1 = 0 et U1 = 0, et à la fin An+1 = 0, Un+1 = U et Ln+1 = L. À chaque étape, les coefficients non triviaux des trois matrices sont stockés dans m : juste avant l’étape k (la k-ème boucle), m(i,j) est égal à
– Ak(i; j) si min(i; j; k) = k,
– Uk(i; j) si min(i; j; k) = i < k,
– Lk(i; j) si min(i; j; k) = j < min(i; k).
Les autres coefficients des trois matrices sont nuls, sauf Lk(i; j) = 1 si i = j < k. Donc (Ak +
LkUk)(i; j) = Ak(i; j) + min(i;j;k 1) Lk(i; l)Uk(l; j).
l=1
Admettons qu’en passant de k à k + 1, on obtient successivement  P
1. Uk+1(i; j) = Uk(i; j) pour tout les (i; j) sauf Uk+1(i; j) = Ak(i; j) pour k = i j,
2. Lk+1(i; j) = Lk(i; j) pour tous les (i; j) sauf Lk+1(i; j) = Ak(i;j) pour k = j i,
Ak(k;k)
3. Ak+1(i; j) = Ak(i; j) [min(i; j; k) = k] Ak(i;k)Ak(k;j) .
Ak(k;k)
On en déduit
min(i;j;k)
Xl
(Ak+1 + Lk+1Uk+1)(i; j) = Ak+1(i; j) + Lk+1(i; l)Uk+1(l; j)
=1
= Ak(i; j) [min(i; j; k) = k]Ak(i; k)Ak(k; j) Ak(k; k)
min(i;j;k 1) Ak(i; k)
+ Lk(i; l)Uk(l; j) + [min(i; j; k) = k] Ak(k; j)
Xl Ak(k; k)
=1
= (Ak + LkUk)(i; j):
C’est la preuve que l’algorithme est correct.
hcalculs 4i
(2)
m(k+1:n,perm(k))=m(k+1:n,perm(k))/m(k,perm(k))
do j=k+1,n
m(k+1:n,perm(j))=m(k+1:n,perm(j))-m(k+1:n,perm(k))*m(k,perm(j))
end do

La résolution LUP x = b

On va résoudre successivement les triangulaires Lz = b et Uy = z, puis P x = y.

Entrées-sorties

On procède à une résolution en place, c’est à dire que le même tableau b contient en entrée le
vecteur b, et en sortie le vecteur x. Les coefficients utiles de L et U sont dans le tableau m, tandis que le tableau trans contient la suite des transpositions de colonnes opérées lors de la décomposition
A = LUP : à l’étape k, les colonnes d’indice k et trans(k) ont été (virtuellement) échangées.
hRésolution LUP x = b 5i
(9) 6.
subroutine lupxeb(m,trans,b)
real, dimension(:,:), intent(in) :: m
integer, dimension(size(m,2)-1), intent(out) :: trans
real, dimension(size(m,1)), intent(inout) :: b
integer, dimension(size(m,2)) :: perm
integer :: n,k,j,ind_piv
real :: r
n=size(m,1)
! m doit etre carree

Initialisation

Il s’agit de restaurer la permutation.
hRésolution LUP x = b 5i+
(9) /5.
perm=(/(k,k=1,n)/)
do k=1,n-1
ind_piv=trans(k) ! enregistre l’indice du k-eme pivot
if (ind_piv/=k) then
! echange perm(k) et perm(ind_piv)
j=perm(k);perm(k)=perm(ind_piv);perm(ind_piv)=j end if
end do

Une descente, une remontée, une permutation

Après la première boucle, b contient le vecteur z. Après la seconde, il contient le vecteur y. Enfin on permute ses éléments pour restaurer x.
7hRésolution LUP x = b 5i+
(9) /6
! resolution Lz=b do k=2,n
b(k:n)=b(k:n)-m(k:n,perm(k-1))*b(k-1) end do
! resolution Uy=z
b(n)=b(n)/m(n,perm(n))
do k=n-1,1,-1
b(1:k)=b(1:k)-m(1:k,perm(k+1))*b(k+1)
b(k)=b(k)/m(k,perm(k))
end do
! resolution de Px=y
do k=n-1,1,-1
ind_piv=trans(k) ! enregistre l’indice du k-eme pivot
if (ind_piv/=k) then ! echange b(k) et b(ind_piv)
r=b(k);b(k)=b(ind_piv);b(ind_piv)=r
end if
end do
end subroutine lupxeb

Produit matrice vecteur

Nous donnons ici un exemple de fonction, qui prend en entrée une matrice et un vecteur, et renvoie en sortie le vecteur résultant de leur produit. Dans une fonction, le résultat est une variable qu’il faut déclarer et qui porte le même nom que la fonction.
Notre fonction pmv est un cas particulier de la fonction prédéfinie matmul. matmul multiplie deux matrices, ou une matrice et un vecteur (ou un vecteur et une matrice). Il est souvent plus efficace de reprogrammer à la main matmul pour le cas particulier que l’on considère.
8 hProduit matrice vecteur 8i (9)
function pmv(a,x)
real, dimension(:,:), intent(in) :: a
real, dimension(size(a,2)), intent(in) :: x
real, dimension(size(a,1)) :: pmv
integer :: n,j
pmv=0.0
n=size(a,2)
do j=1,n
pmv=pmv+a(:,j)*x(j)
end do
end function pmv