Résolution De Système Linéaire - Méthode Directe - JULIA

Le but de ce billet est d’illustrer certaines limites sur le comportement de la résolution de systèmes linéaires à l’aide de méthode directe.

Le langage utilisé est JULIA, l’interface utilisée -le notebook- est Ipython avec le package IJulia. Les graphiques sont réalisés à l’aide du package PyPlot.

Les test seront effectués à l’aide de la méthode fournie dans le package UMFPack (Unsymmetric MultiFrontal) qui est la fonction native de JULIA pour l’opérateur \.

Les différents tests (benchmarks) seront menés sur un système d’équation provenant de la discrétisation différence finie 1d, 2d et 3d du problème modèle :

$$\begin{cases}
-\Delta u = f & \text{dans } \Omega \
u=0& \text{sur } \partial \Omega
\end{cases}$$

Le 15/03/2016

Thierry Clopeau (clopeau@univ-lyon1.fr)


1. Assemblage du système


Laplace 1d

Shéma différence finies sur l’intervalle $[0,1]$ : $\frac{-u_{i+1}+2ui-u{i-1}}{h^2}$ la matrice est tridiagonale.

1ère implémentation

1
2
3
4
5
6
function Laplace1d(n)
h=1./(n+1)
d1=fill(-1/h^2,n-1);
d2=fill(2/h^2,n);
spdiagm((d1,d2,d1),(-1,0,1));
end
Laplace1d (generic function with 1 method)

Laplace 2d

Shéma différence finies sur le domaine $[0,1]^2$ : $\frac{-u{i+n}-u{i+1}+4ui-u{i-1}-u_{i-n}}{h^2}$ la matrice est pentadiagonale on utilise la multiplication de Kronecker

1
2
3
function Laplace2d(n,m)
kron(speye(m),Laplace1d(n))+kron(Laplace1d(m),speye(n))
end
Laplace2d (generic function with 1 method)

Laplace 3d

Shéma différence finies sur le domaine $[0,1]^3$ : $\frac{-u{i+n^2}-u{i+n}-u_{i+1}+6ui-u{i-1}-u{i-n}-u{i-n^2}}{h^2}$ la matrice est heptadiagonale on utilise également la multiplication de Kronecker

1
2
3
function Laplace3d(n,m,p)
kron(speye(p),Laplace2d(n,m))+kron(Laplace1d(p),speye(n*m))
end
Laplace3d (generic function with 1 method)

Test de performance d’assemblage : (matrices de taille max 10 000 000)

Le temps d’assemblage n’est pas forcément optimum, en effet le choix de la fonction kron() n’est pas forcément le meilleur… Néanmoins regardons le temps d’assemblage des matrices qui doit être linéaire en fonction de la taille de la diagonale de la matrice.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
npt=12; ntest=10;
println("1d... ");
N1d=round(Int,logspace(2,7,npt));
T1d=zeros(length(N1d));
Laplace1d(2); #appel pour compiler la fonction
for i=1:length(N1d) tic();for j =1:ntest; Laplace1d(N1d[i]); end;T1d[i]=toq()/ntest ;end;
println(" completed!");

println("2d... ");
N2d=round(Int,logspace(1,3.5,npt));
T2d=zeros(length(N2d));
Laplace2d(2,2); #appel pour compiler la fonction
for i=1:length(N2d) tic();for j =1:ntest; Laplace2d(N2d[i],N2d[i]); end;T2d[i]=toq()/ntest ;end;
println(" completed!");

println("3d...");
N3d=round(Int,logspace(2/3,7/3,npt));
T3d=zeros(length(N3d));
Laplace3d(2,2,2); #appel pour compiler la fonction
for i=1:length(N3d) tic();for j =1:ntest; Laplace3d(N3d[i],N3d[i],N3d[i]); end;T3d[i]=toq()/ntest ;end;
println(" completed!");
1d... 
      completed!
2d... 
      completed!
3d...
      completed!
1
2
3
4
5
6
7
8
using PyPlot
loglog(N1d,T1d,color="red")
loglog(N2d.^2,T2d,color="blue")
loglog(N3d.^3,T3d,color="green")
legend(["1d" ; "2d" ;"3d"],loc=2);
xlabel("taille de la matrice")
ylabel("temps en s")
title("Temps d'assemblage")

png

PyObject <matplotlib.text.Text object at 0x7fe942ad34d0>

En conclusion

A la vue du graphique précédent sur le temps d’assemblage :

  • le 1d (courbe rouge) est visiblement linéaire mais coûte légèrement plus chère que le 2d et 3d !!!!.
  • le 2d (courbe bleue) et le 3d (courbe verte) sont clairement linéaires.

Les courbes obtenues avec surtout l’inversion du 1d laissent à penser que nos méthodes d’assemblage ne sont pas optimales !.

Pour rajouter à cela en terme de performance on pourrait :

  • Utiliser le fait que la matrice soit symétrique pour n’assembler et ne stocker qu’une partie triangulaire.
  • L’usage d’un stockage “sparse” adpaté serait également plus optimal. Par défaut JULIA utilise un format CSC (Compressed Sparse Column), ici un format SKS (Skyline Matrix Storage) serait plus optimal en mémoire et pour le calcul du produit matrice/vecteur.

2. Performance de la bibliothèque de résolution UMFPACK

Dans JULIA (comme MATLAB et d’autres) l’opérateur \ utilise un solveur directe de type décomposition LU multifrontale (voir UMFPACK).

Voici un comparatif du temps de résolution suivant la taille de la matrice :

1
2
3
4
5
6
7
8
9
10
npt=12
println("1d... ");
N1d=round(Int,logspace(2,7,npt));
T1d=zeros(length(N1d));
for i=1:length(N1d)
A=Laplace1d(N1d[i]); b=ones(N1d[i])
println(" n=$(N1d[i])")
tic();A\b;T1d[i]=toq();
end;
println(" completed!");
1d... 
  n=100
  n=285
  n=811
  n=2310
  n=6579
  n=18738
  n=53367
  n=151991
  n=432876
  n=1232847
  n=3511192
  n=10000000
           completed!
1
2
3
4
5
6
7
8
9
println("2d... ");
N2d=round(Int,logspace(1,3.5,npt));
T2d=zeros(length(N2d));
for i=1:length(N2d)
A=Laplace2d(N2d[i],N2d[i]);b=ones(N2d[i].^2)
println(" n=$(N2d[i])")
tic();A\b;T2d[i]=toq();
end;
println(" completed!");
2d... 
  n=10
  n=17
  n=28
  n=48
  n=81
  n=137
  n=231
  n=390
  n=658
  n=1110
  n=1874
  n=3162
          completed!
1
2
3
4
5
6
7
8
9
10
println("3d...");
##### Attention beaucoup de mémoire (8Gb) pour une taille max de 70x70x70 !!!
N3d=round(Int,logspace(2/3,log10(70),npt));
T3d=zeros(length(N3d));
for i=1:length(N3d)
A=Laplace3d(N3d[i],N3d[i],N3d[i]);b=ones(N3d[i].^3)
println(" n=$(N3d[i])")
tic();A\b;T3d[i]=toq();
end;
println(" completed!");
3d...
  n=5
  n=6
  n=8
  n=10
  n=12
  n=16
  n=20
  n=26
  n=33
  n=43
  n=55
  n=70
          completed!
1
2
3
4
5
6
7
8
xlabel("taille de la matrice")
ylabel("temps en s")
title("Temps de résolution UMFPACK")
loglog(N1d,T1d,color="red")
loglog(N2d.^2,T2d,color="blue")
loglog(N3d.^3,T3d,color="green")
legend(["1d" ; "2d" ;"3d"],loc=2)
grid("on")

png

En conclusion

  • On observe un temps d’execution avec des droites en échelle log-log à partir de tailles > $10^4$.
  • En dessous de $10^4$ les courbes sont “flates” et possède un premier point aberrant, cela peut correspondre soit à la précision de la mesure du temps soit à un coût de compilation à la volée (Just In Time) de JULIA.
  • Pour le 1d le coût est linéaire entre $10^4$ et $10^7$.
  • Pour le 2d le coût est légèrement sur-linéaire entre $10^4$ et $5\times10^6$ (théorique $O(N^\frac{3}{2})$.
  • Pour le 3d au delà de $10^4$ le coût se rapproche sérieusement d’un coût quadratique! (théorique $O(N^\frac{5}{3})$

Ici le coût en temps est proportionnel au coût mémoire (non indiqué ici).