]> AND Private Git Repository - ThesisAli.git/blob - Figures/CHAPITRE_02.tex~
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
Update by Ali
[ThesisAli.git] / Figures / CHAPITRE_02.tex~
1 \chapter{R\'esolution de syst\`emes lin\'eaires creux}
2 \label{chap:sys_lin_creux}
3
4 %%-------------------------------------------------------------------------------------------------------%%
5 %Introduction
6 Les systèmes linéaires creux sont utilisés pour la modélisation de nombreux problèmes scientifiques ou industriels,
7 comme la simulation des phénomènes de l'environnement et de la terre ou la transformation industrielle des
8 fluides complexes ou non-newtoniens. De plus, la résolution de ces problèmes fait souvent appel à la résolution
9 de ces systèmes linéaires, qui est considérée comme le processus le plus coûteux en termes de temps de calcul
10 et d'espace mémoire. Par conséquent, la résolution de systèmes linéaires creux doit être aussi efficace que
11 possible, afin de pouvoir traiter des problèmes de taille toujours croissante.
12
13 Il existe, dans le jargon de l'analyse numérique, différentes méthodes de résolution des systèmes linéaires
14 creux. Cependant, le choix d'une méthode de résolution est souvent guidé par les propriétés de la matrice creuse
15 du système à résoudre (structure, symétrie, densité, dimension, etc), la vitesse et la précision
16 de résolution souhaitées. Toutefois, un bon choix de la méthode est aussi guidé par la taille de ces systèmes,
17 qui ralentit davantage le processus de résolution. Par conséquence, ces méthodes doivent être faciles à paralléliser
18 pour exploiter la puissance de calcul et la capacité mémoire des plateformes de calcul parallèle.  
19
20 Dans le présent chapitre, nous décrivons deux classes de méthodes de résolution dans
21 la section~\ref{sec:methodes_resolu}, à savoir les méthodes directes et les méthodes itératives. Dans la
22 section~\ref{sec:stockage}, nous montrons les différentes structures de données des matrices creuses
23 en mémoire, qui permettent d'optimiser l'occupation d'espace mémoire et les accès en lecture/écriture
24 à la mémoire. Enfin dans la section~\ref{sec:parallelisation}, nous présentons la parallélisation des 
25 méthodes itératives pour la résolution de systèmes linéaires creux de grande taille sur les plateformes de 
26 calcul parallèle.
27 %%-------------------------------------------------------------------------------------------------------%%
28
29 %%-------------------------------------------------------------------------------------------------------%%
30 \section{Méthodes de résolution}
31 \label{sec:methodes_resolu}
32 Soit le système linéaire creux suivant:
33 \begin{equation}
34 Ax=b,
35 \label{eq:linear_syst}
36 \end{equation}
37 où $A\in\mathbb{R}^{n\times n}$ est une matrice carrée creuse inversible, $x\in\mathbb{R}^{n}$ est le vecteur
38 solution, $b\in\mathbb{R}^{n}$ est le vecteur du second membre et $n\in\mathbb{N}$ est un grand entier
39 naturel. Les méthodes de résolution des systèmes linéaires sont, plus ou moins, classifiées en deux grandes
40 familles: les méthodes \emph{directes} et les méthodes \emph{itératives}. Dans l'analyse numérique,
41 un algorithme permettant de résoudre un système linéaire creux est appelé un \emph{solveur linéaire creux direct}
42 ou \emph{itératif}, selon le type de la méthode de résolution qu'il utilise.
43
44 \subsection{Méthodes directes}
45 Les méthodes directes sont les méthodes numériques les plus anciennes pour la résolution des systèmes
46 linéaires. De plus, jusqu'aux années quatre-vingts,
47 %dans les années précédentes, 
48 elles ont souvent été préférées aux méthodes itératives,
49 en raison de leur \emph{robustesse} et de leur comportement \emph{prévisible}. En effet, que ce soit
50 pour des matrices denses ou creuses, ces méthodes permettent d'obtenir, théoriquement en l'absence des erreurs
51 d'arrondi, une solution \emph{exacte} en un nombre d'opérations élémentaires \textit{fini} et \emph{connu a priori}.
52 La plupart des méthodes directes sont basées, principalement, sur deux étapes de calcul: la factorisation
53 de la matrice suivant ses propriétés (symétrie, définie positive, etc) et la résolution
54 par des substitutions successives~\cite{ref1}. 
55
56 Parmi les solveurs directs les plus utilisés, nous pouvons citer le \textit{solveur LU} basé sur la
57 méthode de factorisation \textit{LU}~\cite{ref2}. Dans un premier temps, ce solveur procède à la factorisation de
58 la matrice des coefficients $A$ du système linéaire à résoudre (\ref{eq:linear_syst}). Dans cette étape,
59 la matrice $A$ est décomposée en deux matrices triangulaires inférieure $L=(l_{ij})_{1\leqslant i,j\leqslant n}$
60 et supérieure $U=(u_{ij})_{1\leqslant i,j\leqslant n}$, de façon à ce que la matrice $A$ puisse s'écrire sous
61 forme d'un produit des deux matrices $L$ et $U$. Une fois la factorisation établie, le solveur \textit{LU}
62 résout le système linéaire en deux étapes:
63 \begin{itemize*}
64 \item la résolution du système $Ly=b$ par l'algorithme de \emph{descente} (substitutions successives
65 avant ou résolution de haut vers le bas),
66 \item puis, la résolution du système $Ux=y$ par l'algorithme de \emph{remontée} (substitutions successives
67 arrière ou résolution de bas vers le haut).
68 \end{itemize*}
69
70 Cependant, les méthodes directes de résolution des systèmes linéaires creux sont souvent plus
71 compliquées que celles des systèmes denses. La principale difficulté réside dans l'étape de factorisation
72 de la matrice. En effet, durant la construction des deux facteurs triangulaires $L$ et $U$, certains
73 zéros de la matrice creuse d'origine, $A$, peuvent devenir des valeurs non nulles dans les matrices $L$
74 et $U$, ce que l'on appelle, communément, processus de \emph{remplissage}. Par conséquent, une mise en \oe uvre
75 efficace d'un solveur direct creux dépend largement des structures de données en mémoire (formats de stockage)
76 des valeurs non nulles des matrices $L$ et $U$. De ce fait, un solveur direct creux est généralement
77 basé sur quatre étapes distinctes~\cite{ref3,ref4} au lieu de deux comme dans les solveurs denses:
78 \begin{enumerate*}
79  \item Réorganisation de la matrice $A$: qui consiste à réordonner les lignes et les colonnes de $A$,
80 et donc les valeurs non nulles, de telle sorte que les facteurs $L$ et $U$ subissent peu de remplissage
81 à l'étape de factorisation numérique;
82  \item Factorisation symbolique: qui détermine a priori les structures de données appropriées aux valeurs
83 non nulles des facteurs $L$ et $U$;
84  \item Factorisation numérique: qui calcule les facteurs $L$ et $U$;
85  \item Résolution du système linéaire creux: qui effectue les substitutions avant-arrière en utilisant
86 les deux facteurs $L$ et $U$ calculés à l'étape précédente. 
87 \end{enumerate*}
88 Dans le cas des matrices creuses symétriques et définies positives (valeurs propres strictement posistives),
89 les quatres étapes peuvent être bien
90 séparées comme dans le solveur \textit{CHOLMOD}~\cite{ref5}. Pour les matrices creuses asymétriques,
91 les solveurs directs peuvent fusionner les étapes 2 et 3 en une seule étape comme les solveurs \textit{SuperLU}~\cite{ref6}
92 et \textit{MUMPS}~\cite{ref20},
93 ou même fusionner les étapes 1, 2 et 3 comme le solveur \textit{UMFPACK}~\cite{ref7}.
94
95 Toutefois, les méthodes directes ne sont pas bien adaptées pour la résolution des systèmes linéaires
96 creux de très grandes dimensions, et en particulier, ceux associés aux problèmes tridimensionnels. En fait,
97 l'occupation espace mémoire du processus de factorisation croît \textit{superlinéairement}
98 avec la taille de la matrice et dépend fortement de la structure de cette dernière (positions des valeurs non nulles)
99 pour réduire le remplissage. Hélas, quelque soit le nombre de fois que les coefficients de la
100 matrice sont réordonnés, cela n'empêche pas le remplissage de croître rapidement et proportionnellement
101 à la taille de la matrice. En outre, les solveurs directs sont très coûteux en termes de temps de résolution
102 qui est de $\mathcal{O}(n^{3})$, où $n$ représente la taille du système. Dans ce cas, les méthodes itératives
103 décrites dans la section~\ref{sec:iterative}, sont particulièrement plus utiles lorsque les systèmes linéaires
104 à résoudre sont creux et de très grande dimension car, contrairement aux méthodes directes, aucun processus de remplissage n'est nécessaire.  
105
106 \subsection{Méthodes itératives}
107 \label{sec:iterative}
108 Les méthodes itératives procédent par itérations successives
109 d'un même bloc d'opérations élémentaires, au cours desquelles une séquence infinie de solutions approximatives
110 $\{x^{(k)}\}_{k\geq0}$ est générée. A partir d'une solution initiale $x^{(0)}$ (choisie pour la première
111 itération), une méthode itérative détermine à chaque itération
112 $k\geq1$, une solution approximative, $x^{(k)}$, qui converge progressivement vers la solution exacte, $x^{*}$,
113 du système linéaire~\ref{eq:linear_syst}, telle que:
114 \begin{equation}
115 x^{*}=\lim\limits_{k\to\infty}x^{(k)}=A^{-1}b\in\mathbb{R}^{n}.
116 \end{equation}
117
118 A la différence d'une méthode directe, le nombre d'itérations d'une méthode itérative nécessaire pour
119 atteindre la solution exacte du système linéaire à résoudre \textit{n'est pas connu} au préalable et
120 peut être \textit{infini}. Cependant, dans la pratique, un solveur itératif donne souvent une solution
121 approximative satisfaisant une précision requise et/ou après un nombre maximum d'itérations fixé a priori.
122 Les critères d'arrêt les plus fréquemment utilisés par les solveurs itératifs sont les suivants:
123 \begin{itemize}
124   \item $\|x^{(k)}-x^{(k-1)}\|<\varepsilon$: erreur absolue,
125   \item $\|x^{(k)}-x^{(k-1)}\|<\varepsilon\|x^{(k-1)}\|$: erreur relative,
126   \item $\|b-Ax^{(k)}\|<\varepsilon\|b\|$: résidu de la solution. 
127 \end{itemize}
128 où $x^{(k)}$ est la solution approximative obtenue à l'itération $k\leqslant K$, $K$ étant le nombre maximum
129 d'itérations, et $\varepsilon$ est le seuil de tolérance requis pour la précision de calcul. Ainsi, un solveur
130 linéaire itératif est dit \textit{convergent} si pour tout choix de la solution initiale $x^{(0)}\in\mathbb{R}^{n}$,
131 l'erreur ou le résidu de la solution finale $x^{(k)}\in\mathbb{R}^{n}$, obtenue après $k$ itérations, est
132 suffisamment \textit{petit}. Cependant, la convergence des méthodes itératives n'est pas toujours assurée
133 et dépend fortement des caractéristiques du système linéaire, à savoir: la structure creuse, le conditionnement
134 et le rayon spectral de la matrice des coefficients $A$. Il existe différentes classes de méthodes itératives
135 et le choix de l'utilisation d'une méthode est souvent guidé par son efficacité de résolution et les propriétés du système
136 linéaire à résoudre. Nous citons ici, globalement, trois classes de méthodes itératives:
137 stationnaires, non stationnaires et multigrilles.
138
139 \subsubsection{Méthodes stationnaires}
140 Les méthodes itératives \textit{stationnaires} permettant de résoudre le système linéaire (\ref{eq:linear_syst}),
141 peuvent être exprimées sous la forme simple suivante:
142 \begin{equation}
143 \forall k\in\mathbb{N}, x^{(k+1)}=f(x^{(k)})=Cx^{(k)}+d,
144 \end{equation}
145 telle que $f$ soit une fonction \textit{affine} représentant un schéma de point fixe, où la matrice d'itération $C$
146 et le vecteur $d$ restent \textit{inchangés} pour toutes les itérations $k$. Les méthodes itératives
147 stationnaires sont basées sur la décomposition de la matrice des coefficients, $A$, comme suit:
148 \begin{equation}
149 A=M-N,
150 \end{equation}
151 de sorte que le système linéaire (\ref{eq:linear_syst}) puisse s'écrire sous forme: $Mx=Nx+b$ à laquelle est 
152 associée l'itération suivante:
153 \begin{equation}
154 x^{(k+1)}=Cx^{(k)}+d=M^{-1}Nx^{(k)}+M^{-1}b,
155 \end{equation}
156 où $M$ est une matrice inversible. Par conséquent, une méthode itérative stationnaire est \textit{convergente}
157 vers la solution exacte $x^{*}$ si, et seulement si, le rayon spectral (la plus grande valeur propre) de la matrice $C$ est strictement
158 inférieur à $1$: $\rho(C)<1$. Les principales méthodes itératives stationnaires sont basées sur la décomposition
159 de la matrice des coefficients en matrices diagonale $D$, triangulaire strictement inférieure $L$ et
160 triangulaire strictement supérieure $U$: $A=D-L-U$. Parmi ces méthodes itératives nous pouvons citer~\cite{ref8}:
161 \begin{itemize}
162   \item Jacobi: $M=D$ et $N=L+U$,
163   \item Gauss-Seidel: $M=D-L$ et $N=U$, 
164   \item Surrelaxation successive (SOR): $M=D-\omega L$ et $N=\omega U+(1-\omega)D$ et $\omega>0$, 
165   \item Richardson: $C=I-\omega A$, $d=\omega b$ et $\omega>0$.
166 \end{itemize}
167
168 L'algorithme~\ref{alg:stationnaire} présente les principaux points clés d'un solveur itératif
169 stationnaire. Il permet de calculer, à partir d'une solution initiale $x^{(0)}$, une série de
170 solutions approximatives successives $x$, jusqu'à ce que la convergence soit atteinte
171 (dans l'algorithme, à la ligne $7$, résidu<$\varepsilon$).
172
173 \begin{algorithm}[H]
174   \SetLine
175   \linesnumbered
176   \Entree{$A$ (matrice), $C$ (matrice d'itération), $b$ (vecteur), $d$ (vecteur)}
177   \Sortie{$x$ (vecteur)}
178      
179   \BlankLine
180
181   choisir $x^{(0)}$ solution initiale, $\varepsilon$ seuil de tolérance\;
182   variables $k$, $convergence$\;
183   $convergence \leftarrow faux$\;
184   $k \leftarrow 0$\;
185   
186   \Tq{$\neg convergence$}
187   {
188     évaluer: $x^{(k+1)}=f(x^{(k)})=Cx^{(k)}+d$\;
189     $convergence \leftarrow (\|b-Ax^{(k+1)}\|\cdot\|b\|^{-1}<\varepsilon)$\;
190     $k \leftarrow k+1$\;
191   }
192 \caption{Algorithme général d'un solveur itératif stationnaire}
193 \label{alg:stationnaire}
194 \end{algorithm}
195
196 La ligne $6$ de l'algorithme~\ref{alg:stationnaire} permet de calculer les composantes du
197 vecteur $x$ en appliquant la fonction affine de la méthode itérative stationnaire utilisée.
198 Par exemple, dans les méthodes Jacobi et Gauss-Seidel, les $n$ composantes du vecteur
199 $x$ à l'itération $(k+1)$ sont calculées comme suit: 
200
201 \begin{tabular}{rl}
202   $Jacobi:$& $ x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$\hspace{1.6cm}$\Leftrightarrow$\\ \\
203   & $x^{(k+1)}_{i}=(b_{i}-\displaystyle\sum_{j\neq i}a_{ij}x^{(k)}_{i})/a_{ii}$, $\forall i,j$, $ 1\leqslant i,j\leqslant n$,  \\ \\
204   $Gauss-Seidel:$ & $x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b$\hspace{0.5cm}$\Leftrightarrow$\\ \\
205   & $x^{(k+1)}_{i}=(b_{i}-\displaystyle\sum_{j<i}a_{ij}x^{(k+1)}_{i}-\displaystyle\sum_{j>i}a_{ij}x^{(k)}_{i})/a_{ii}$, $\forall i,j$, $ 1\leqslant i,j\leqslant n$. 
206 \end{tabular}
207
208 Comme nous l'avons évoqué, précédemment, le principal critère de choix entre les solveurs
209 itératifs est l'efficacité, et donc la convergence, des méthodes de résolution utilisées
210 par ces derniers. Par exemple, pour trouver la nouvelle solution $x$ à l'itération $(k+1)$,
211 la méthode Jacobi utilise toutes les composantes du vecteur $x$ calculées à l'itération
212 précédente $k$, alors que la méthode Gauss-Seidel utilise les nouvelles composantes
213 calculées à l'itération courante $(k+1)$ dès que possible. Par conséquent, cette technique
214 de mise à jour des composantes du vecteur solution $x$ permet au solveur Gauss-Seidel
215 d'effectuer moins d'itérations et donc d'avoir, en général, une convergence plus rapide que
216 le solveur Jacobi.
217
218 Les méthodes stationnaires sont les méthodes itératives les plus anciennes et leurs
219 algorithmes sont simples à comprendre et à mettre en \oe uvre. Cependant, elles sont souvent
220 moins efficaces pour la résolution de beaucoup de systèmes linéaires creux, vu leur applicabilité
221 limitée à certains types de matrices. Ces dernières années, les méthodes non stationnaires,
222 définies dans la section suivante, sont devenues les méthodes itératives les plus utilisées
223 pour leur efficacité à résoudre plusieurs types de systèmes linéaires.
224
225 \subsubsection{Méthodes non stationnaires}
226 \label{sec:non-stationnaire}
227 Les méthodes \textit{non stationnaires} sont les méthodes itératives les 
228 plus récentes; leurs algorithmes sont généralement plus difficiles à comprendre et à 
229 mettre en \oe uvre, mais elles peuvent être plus efficaces pour la résolution de certains systèmes
230 linéaires creux. Elles peuvent être exprimées sous la forme suivante:
231 \begin{equation}
232   \forall k\in\mathbb{N}, x^{(k+1)}=f_{k}(x^{(k)}),
233 \end{equation}
234 où, contrairement aux méthodes stationnaires, la fonction $f_{k}$ permettant de calculer
235 le nouvel itéré $x^{(k+1)}$ diffère d'une itération à l'autre. En fait, les méthodes non stationnaires 
236 impliquent des données qui \textit{changent} à chaque itération et qui déclinent des opérations
237 vectorielles utilisées par les méthodes itératives, comme les produits scalaires et les 
238 opérations \textit{AXPY} ($y\leftarrow ax+y$).
239
240 Le principe général des méthodes non stationnaires est la construction d'une séquence de 
241 vecteurs orthogonaux et les projections sur des sous-espaces. L'une des classes importantes
242 des méthodes non stationnaires est celle des \textit{méthodes itératives de Krylov}.
243 L'idée générale d'une méthode de Krylov~\cite{ref8} est de chercher une approximation, $x^{(k)}$, de 
244 la solution de système linéaire (\ref{eq:linear_syst}) dans un sous-espace de Krylov,
245 $\mathcal{K}_{k}$, d'ordre $k$:
246 \begin{equation}
247 x^{(k)}-x^{(0)}=q_{k-1}(A)r^{(0)}\in\mathcal{K}_{k}(A,r^{(0)}),
248 \end{equation}
249 telle que la condition Petrov-Galerkin soit satisfaite:\[b-Ax^{(k)}\bot\mathcal{L}_{k},\]
250 où $q_{k-1}$ est un polynôme de degré $k-1$ et $\mathcal{L}_{k}$ est un autre sous-espace
251 vectoriel d'ordre $k$. Le sous-espace de Krylov d'ordre $k$, noté $\mathcal{K}_{k}(A,r^{(0)})$, 
252 est un sous-espace vectoriel construit sur $A$ et $r^{(0)}$ et engendré par les $k$ vecteurs
253 orthogonaux deux à deux $r^{(0)}$, $Ar^{(0)}$, $\dots$, $A^{(k-1)}r^{(0)}$ (base orthonormée):
254 \[\mathcal{K}_{k}(A,r^{(0)})\equiv \text{engendré}\{r^{(0)}, Ar^{(0)}, A^{2}r^{(0)}, \dots, A^{k-1}r^{(0)}\},\]
255 où $r^{(0)}=b-Ax^{(0)}$ est le résidu initial de la solution initiale $x^{(0)}$ choisie, plus 
256 ou moins aléatoirement, au démarrage du solveur de Krylov. 
257
258 La plupart des méthodes de Krylov sont basées sur des itérations ne requierant que des multiplications
259 de la matrice du système par un vecteur, des produits scalaires et des additions de vecteurs; et
260 elles diffèrent généralement dans le choix de $\mathcal{L}_{k}$, le calcul d'une base orthonormée
261 pour $\mathcal{K}_{k}$ (certaines méthodes utilisent la matrice creuse $A$ et d'autres utilisent 
262 sa transposée $A^{T}$) et la réinitialisation de cette base à toutes les $m$ itérations (restarts).
263 Les méthodes de Krylov ont, généralement, une convergence plus rapide que les méthodes stationnaires.
264 Néanmoins, pour certains types de matrices creuses, elles n'arrivent pas à converger ou elles
265 effectuent de nombreuses itérations pour converger. En pratique, elles sont, généralement,
266 utilisées avec un \textit{préconditionneur} qui permet d'améliorer et/ou d'accélérer leur convergence.
267 Le processus de préconditionnement consiste à remplacer le système linéaire creux (\ref{eq:linear_syst})
268 par l'un des deux systèmes suivants:
269 \begin{equation}M^{-1}Ax=M^{-1}b,\label{eq:left}\end{equation}
270 \begin{center}ou\\\end{center}
271 \begin{equation}\begin{array}{cc}AM^{-1}\hat x=b, & x=M^{-1}\hat x,\end{array}\label{eq:right}\end{equation}
272 qui sont plus faciles à résoudre. $M$ représente la matrice de préconditionnement et les formules (\ref{eq:left})
273 et (\ref{eq:right}) représentent, respectivement, le préconditionnement à gauche et le préconditionnement
274 à droite du système linéaire $Ax=b$. Il existe beaucoup de méthodes de Krylov~\cite{ref8,ref9}, par exemple
275 \textit{CG} ($\mathcal{L}_{k}=\mathcal{K}_{k}$), \textit{GMRES} ($\mathcal{L}_{k}=A\mathcal{K}_{k}$) et \textit{BICG}
276 ($\mathcal{L}_{k}=\mathcal{K}_{k}(A^{T},r^{(0)})$).
277
278 \subsubsection{Méthodes multigrilles}
279 Une autre catégorie de méthodes itératives est celle des méthodes \textit{multigrilles} dont les origines
280 remontent aux années soixante \cite{ref10}. Initialement conçues pour la résolution de l'équation de Poisson,
281 les méthodes multigrilles ont rapidement évolué pour être appliquées à un large éventail d'applications dans
282 le domaine de la physique ou des mathématiques des finances, comme par exemple la résolution des problèmes
283 linéaires ou non linéaires avec des conditions aux limites \cite{ref11} ou des problèmes d'optimisation \cite{ref12}.
284
285 Les méthodes itératives standards, comme Jacobi et Gauss-Seidel, sont caractérisées par un faible taux de
286 convergence global, $\tau$, qui peut être défini à l'itération $k$ comme suit:
287 \begin{equation}
288   \tau_{k}=\frac {\|r^{(k)}\|} {\|r^{(k-1)}\|} = \frac {\|b-Ax^{(k)}\|} {\|b-Ax^{(k-1)}\|},
289 \end{equation}
290 où $r^{(k)}$ est le résidu de la solution approximative $x^{(k)}$ trouvée après $k$ itérations. En fait, malgré
291 leur convergence rapide aux premières itérations, ces méthodes itératives peinent à maintenir leur vitesse
292 de convergence tout le long de la résolution. Cela peut être expliqué par leur incapacité à réduire les basses
293 fréquences de l'erreur, ce qui se traduit par la décroissance lente de la norme du résidu durant les dernières
294 itérations de la résolution. Par conséquent, beaucoup d'itérations sont nécessaires pour donner une solution
295 approchée intéressante lorsque la vitesse de convergence est faible. A cet effet, les méthodes multigrilles
296 sont conçues pour améliorer la vitesse de convergence de ces méthodes itératives et ainsi, réduire leur coût
297 de résolution. 
298
299 Comme leur nom l'indique, les méthodes multigrilles utilisent différents niveaux de maillages réguliers
300 (au moins deux grilles), du plus fin au plus grossier, chaque n\oe ud de maillage étant espacé d'un autre
301 d'une distance $2qh$, tel que $q\in\mathbb{N}$ représente le niveaux de la grille. Pour des raisons de
302 simplicité, nous ne décrirons, ci-après, que le principe général de résolution d'une méthode à deux grilles
303 sur laquelle sont basées, d'ailleurs, l'ensemble des méthodes multigrilles. On considère un système
304 linéaire creux (\ref{eq:fine}) discrétisé sur un maillage $\Omega_{h}$ de pas $h=\frac{1}{n+1}$, $n$
305 étant le nombre de points de discrétisation:
306 \begin{equation}
307   \begin{array}{lll}
308     A_{h}x_{h}=b_{h}, & sur & \Omega_{h},\\
309   \end{array} 
310 \label{eq:fine}
311 \end{equation}
312 où $A_{h}$, $x_{h}$ et $b_{h}$ représentent, respectivement, la matrice creuse, le vecteur solution et le
313 second membre du système discrétisé sur $\Omega_{h}$. La résolution de ce système par une méthode itérative
314 permet de trouver seulement une approximation $x_{h}$ à la solution exacte $x_{h}^{*}$ sur $\Omega_{h}$.
315 L'erreur algébrique $e_{h}$ et le résidu $r_{h}$ de cette solution approximative peuvent être calculés comme
316 suit: 
317 \begin{equation}
318   e_{h}=x^{*}_{h}-x_{h},
319 \label{eq:erreur}
320 \end{equation}
321 \begin{equation}
322   r_{h}=b_{h}-A_{h}x_{h},
323 \label{eq:residu}
324 \end{equation}
325 De ces deux dernières formules, nous pouvons déduire l'équation de l'erreur suivante:
326 \begin{equation}
327  \begin{array}{lll}
328   A_{h}e_{h}=r_{h}, & sur & \Omega_{h},
329  \end{array} 
330 \label{eq:residuelle}
331 \end{equation} 
332 L'objectif de la méthode bigrille \cite{ref13} est de réduire les basses fréquences de l'erreur sur
333 la grille grossière $\Omega_{2h}$, que la méthode itérative n'a pas réussie à éliminer sur la gille fine
334 $\Omega_{h}$. En fait, les basses fréquences de l'erreur sur la grille fine apparaissent comme des hautes fréquences de
335 l'erreur sur la grille grossière, qui sont faciles à éliminer par les méthodes itératives. Le principe
336 général d'une méthode bigrille repose sur la résolution de l'équation de l'erreur (\ref{eq:residuelle}) sur la grille
337 grossière $\Omega_{2h}$, pour trouver une erreur $e_{h}$ permettant de corriger la valeur de la solution
338 approximative $x_{h}$: $x_{h}=x_{h}+e_{h}$, obtenue après quelques itérations sur la grille fine
339 $\Omega_{h}$. L'algorithme \ref{alg:bigrille} montre les principaux points clés d'une méthode bigrille.
340 Les processus de pré-lissage et de post-lissage consistent à effectuer quelques itérations d'un solveur
341 itératif (par exemple Jacobi ou Gauss-Seidel) sur la grille fine, afin d'éliminer les hautes fréquences
342 de l'erreur. Les deux opérateurs de restriction $R^{2h}_{h}$ et de prolongement $P^{h}_{2h}$ permettent le
343 transfert des données entre la grille fine et la grille grossière. $A_{2h}$ est la matrice d'itération sur
344 la grille grossière, telle que: $A_{2h}=R_{h}^{2h}A_{h}P_{2h}^{h}$. La résolution de l'équation de l'erreur
345 avec une méthode itérative sur la grille grossière (à la ligne $7$ de l'algorithme \ref{alg:bigrille})
346 permet de trouver l'erreur de correction de la solution approximative $x_{h}$ obtenue sur la grille fine
347 et ainsi, d'éliminer les basses fréquences de l'erreur. 
348
349 \begin{algorithm}[H]
350   \SetLine
351   \linesnumbered
352   \Entree{$A$ (matrice), $b$ (vecteur)}
353   \Sortie{$x$ (vecteur)}
354      
355   \BlankLine
356
357   choisir $x_{h}^{0}$ solution initiale, $\varepsilon$ seuil de tolérance\;
358   $convergence \leftarrow faux$\;
359   
360   \Tq{$\neg convergence$}
361   {
362     Pré-lissage: résoudre $A_{h}x_{h}=b_{h}$ sur $\Omega_{h}$\;
363     Calcul du résidu: $r_{h}=b_{h}-A_{h}x_{h}$\;
364     Restriction sur $\Omega_{2h}$: $r_{2h}=R_{h}^{2h}r_{h}$\;
365     Résoudre: $A_{2h}e_{2h}=r_{2h}$ sur $\Omega_{2h}$\;
366     Prolongement sur $\Omega_{h}$: $e_{h}=P^{h}_{2h}e_{2h}$\;
367     Correction: $x_{h}=x_{h}+e_{h}$\;
368     Post-lissage: résoudre $A_{h}x_{h}=b_{h}$ sur $\Omega_{h}$\; 
369     $convergence \leftarrow (\|b_{h}-A_{h}x_{h}\|<\varepsilon)$\;
370     $x_{h}^{0}\leftarrow x_{h}$\;
371   }
372 \caption{Algorithme général d'une méthode bigrille}
373 \label{alg:bigrille}
374 \end{algorithm}
375
376 Cependant, dans la pratique, les méthodes multigrilles ne se limitent pas à deux grilles
377 pour obtenir une convergence rapide. Elles sont, souvent, appliquées avec plusieurs grilles
378 de plus en plus grossières: $\Omega_{h}$, $\Omega_{2h}$, $\Omega_{4h}$, $\Omega_{8h}$, $\cdots$.
379 L'idée générale de ces méthodes est l'application récursive du principe de la méthode bigrille.
380 En effet, les méthodes multigrilles appliquent une méthode bigrille pour résoudre l'équation
381 de l'erreur (ligne 7 de l'algorithme~\ref{alg:bigrille}), afin de trouver une meilleure correction
382 de la solution approximative. Il existe plusieurs versions de ce type de méthodes, par exemple
383 \textit{V-cycle}, \textit{W-cycle} ou \textit{FMG}, qui sont présentées dans~\cite{ref14,ref15,ref16}.
384 %%-------------------------------------------------------------------------------------------------------%%
385
386 %%-------------------------------------------------------------------------------------------------------%%
387 \section{Formats de stockage des matrices creuses}
388 \label{sec:stockage}
389 La plupart des méthodes itératives sont basées sur des opérations de calcul traitant des vecteurs et des
390 matrices (voir section~\ref{sec:iterative}). Cependant, parmi toutes ces opérations, seule la multiplication
391 matrice-vecteur est considérée comme l'opération la plus importante des méthodes itératives, car elle est
392 souvent très coûteuse en termes de temps d'exécution et d'occupation espace mémoire. Par conséquent,
393 l'efficacité d'une méthode itérative dépend fortement de la performance de calcul de sa multiplication
394 matrice-vecteur.
395
396 En outre, dans le cas des systèmes linéaires creux, la multiplication matrice-vecteur nécessite, aussi, une
397 attention très particulière pour le format de stockage de la matrice creuse dans la mémoire. Le stockage
398 naïf, ligne par ligne ou colonne par colonne, d'une matrice creuse peut causer une perte considérable
399 d'espace mémoire et de temps de calcul. En fait, les éléments nuls de la matrice sont sauvegardés dans
400 la mémoire alors que leur prise en compte dans la multiplication est totalement inutile. De plus, la
401 nature creuse de la matrice engendre, souvent, des accès irréguliers à la mémoire pour la lecture de ses
402 éléments non nuls et ceci, ne peut que dégrader encore davantage la performance de la multiplication
403 matrice-vecteur.
404
405 Néanmoins, une matrice creuse permet, aussi, de tirer profit de son aspect creux et ce, en ne sauvegardant
406 que ses éléments non nuls, qui permettra ainsi, de réduire la taille de la mémoire nécessaire
407 pour le stockage et d'accélérer le calcul de la multiplication matrice-vecteur. Il existe dans la
408 littérature~\cite{ref8,ref17}, plusieurs formats de stockage compressés pour les matrices creuses. Ils permettent un accès
409 facile aux éléments non nuls en connaissant leurs coordonnées dans la mémoire et une exécution intelligente
410 des opérations de calcul sur ces éléments. Dans la présente section, nous décrivons au total cinq
411 formats de compression des matrices creuses, qui sont intéressants pour la suite de ce document.
412
413 \subsection{COO}
414 \textit{COO} (Coordinate) est un format de stockage particulièrement simple. Il permet de représenter une matrice
415 creuse sous forme d'une structure de données composée de trois tableaux: un tableau de réels, $Val$, contenant
416 tous les éléments non nuls de la matrice creuse $A$ et deux tableaux d'entiers,
417 $Lig$ et $Col$, contenant les indices de lignes et les indices de colonnes de ces éléments non nuls, respectivement.
418 Les trois tableaux ont une taille $nnz$ (nombre d'éléments non nuls de la matrice). La figure~\ref{fig:coo}
419 décrit le format \textit{COO} pour une matrice creuse d'ordre $(4\times4)$, $A$, et l'algorithme~\ref{alg:coo} présente
420 la multiplication matrice-vecteur pour ce format de stockage. 
421 \begin{figure}[!h]
422 \begin{center}
423 \begin{tabular}{c}
424 $A=\left[
425 \begin{matrix}
426 a_{00} & 0 & 0 & a_{03} \\
427 0 & 0 & a_{12} & 0 \\
428 a_{20} & a_{21} & a_{22} & 0 \\
429 0 & a_{31} & 0 & a_{33} \\
430 \end{matrix}
431 \right]
432 $
433 \\
434 \\
435 $
436 \begin{array}{cccccccccccc}
437 Val & = & [ & a_{00} & a_{03} & a_{12} & a_{20} & a_{21} & a_{22} & a_{31} & a_{33} & ] \\
438 Lig & = & [ & 0     & 0      & 1     & 2      & 2     & 2      & 3     & 3      & ] \\
439 Col & = & [ & 0     & 3      & 2     & 0      & 1     & 2      & 1     & 3      & ] \\
440 \end{array}
441 $
442 \end{tabular}
443 \caption{Le format de stockage COO}
444 \label{fig:coo}
445 \end{center}
446 \end{figure}
447
448 Dans le format \textit{COO}, les éléments non nuls d'une matrice creuse peuvent être stockés ligne par ligne (voir
449 l'exemple de la figure~\ref{fig:coo}), colonne par colonne ou dans un ordre quelconque. Cependant, ce format
450 n'est efficace que pour les matrices creuses non structurées, autrement, nous pouvons remarquer que dans le
451 cas de stockage ligne par ligne (ou colonne par colonne), le tableau $Lig$ (respectivement, le tableau $Col$)
452 peut contenir des informations redondantes successives qui pointent le même indice de ligne (respectivement, le même indice de
453 colonne). Cela impliquerait un espace mémoire non négligeable pour stocker une grande matrice creuse structurée.
454
455 \begin{algorithm}
456   \SetLine
457   \linesnumbered
458   \Entree{$Val$, $Lig$ et $Col$ (matrice), $nnz$ (nombre d'éléments non nuls), $n$ (taille de la matrice), $x$ (vecteur)}
459   \Sortie{$y$ (vecteur)}
460   \BlankLine
461   variables $i$, $j$, $id$\; 
462   \Pour{$i=0$ \KwA $n-1$}{
463     $y[i] \leftarrow 0$\;
464   }
465   \BlankLine
466   \Pour{$id=0$ \KwA $nnz-1$}{
467     $i \leftarrow Lig[id]$\;
468     $j \leftarrow Col[id]$\;
469     $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
470   }
471 \caption{Multiplication matrice-vecteur avec le format COO}
472 \label{alg:coo}
473 \end{algorithm}
474
475 \subsection{CSR vs. CSC}
476 Les formats \textit{CSR} (Compressed Sparse Row) et \textit{CSC} (Compressed Sparse Column) sont, probablement,
477 les plus connus et les plus utilisés pour le stockage des matrices creuses. En effet, ce
478 sont deux formats qui ne posent aucune condition sur la structure creuse de la matrice et
479 ne permettent de stocker que les informations nécessaires sur la matrice.
480 Ils permettent de représenter la matrice creuse sous forme de trois tableaux de données. Comme
481 le format \textit{COO}, le format \textit{CSR} (respectivement, \textit{CSC}) stocke les $nnz$ éléments non nuls de la matrice
482 et leurs indices de colonnes (respectivement, leurs indices de lignes) ligne par ligne (respectivement, colonne
483 par colonne) dans deux tableaux distincts $Val$ et $Col$ (respectivement, $Val$ et $Lig$). Un troisième
484 tableau d'une taille $n+1$, $Ptr$, est utilisé pour stocker la position du premier élément
485 non nul de chaque ligne de la matrice (respectivement, de chaque colonne de la matrice). La figure~\ref{fig:csr}
486 donne les formats de stockage \textit{CSR} et \textit{CSC} de la matrice creuse $A$ donnée comme exemple.
487 Les algorithmes~\ref{alg:csr} et~\ref{alg:csc} présentent une multiplication
488 matrice-vecteur avec le format \textit{CSR} et celle avec le format \textit{CSC}, respectivement.
489
490 \begin{figure}[!h]
491 \begin{center}
492 \begin{tabular}{c}
493 $A=\left[
494 \begin{matrix}
495 a_{00} & 0      & 0     & a_{03} \\
496 0     & 0      & a_{12} & 0 \\
497 a_{20} & a_{21} & a_{22} & 0 \\
498 0     & a_{31} & 0      & a_{33} \\
499 \end{matrix}
500 \right]
501 $
502 \\
503 \\
504 $
505 CSR:\left\{
506 \begin{array}{ccccccccccccc}
507 Val & = & [ & a_{00} & a_{03} & a_{12} & a_{20} & a_{21} & a_{22} & a_{31} & a_{33} &   & ] \\
508 Col & = & [ & 0     & 3      & 2     & 0      & 1     & 2      & 1     & 3      &   & ] \\
509 Ptr & = & [ & 0     &        & 2     & 3      &       &        & 6     &        & 8 & ] \\ 
510 \end{array}
511 \right.
512 $
513 \\
514 \\
515 $
516 CSC:\left\{
517 \begin{array}{ccccccccccccc}
518 Val & = & [ & a_{00} & a_{20} & a_{21} & a_{31} & a_{12} & a_{22} & a_{03} & a_{33} &   & ] \\
519 Lig & = & [ & 0     & 2      & 2     & 3      & 1     & 2      & 0     & 3      &   & ] \\
520 Ptr & = & [ & 0     &        & 2     &        & 4     &        & 6     &        & 8 & ] \\ 
521 \end{array}
522 \right.
523 $
524 \end{tabular}
525 \caption{Le format de stockage CSR et CSC}
526 \label{fig:csr}
527 \end{center}
528 \end{figure}
529
530 Les formats \textit{CSR} et \textit{CSC} sont sûrement les plus économes en espace mémoire pour
531 le stockage des matrices creuses mais, ils ne sont pas forcément les plus efficaces.
532 En effet, d'après l'algorithme~\ref{alg:csr}, chaque opération atomique ($y[i] \leftarrow y[i] + Val[id] \times x[j]$) de la
533 multiplication matrice-vecteur \textit{CSR} nécessite le chargement de l'indice de colonne
534 et un accès indirect à la mémoire pour lire la valeur du vecteur $x$ correspondante.
535 De même, pour le format \textit{CSC} (algorithme~\ref{alg:csc}), chaque opération atomique
536 de la multiplication requiert un chargement de l'indice de ligne et puis un accès
537 indirect à la mémoire pour écrire la valeur du vecteur $y$ à l'adresse
538 correspondante. Donc, la multiplication sur des matrices ayant des lignes ou des
539 colonnes très creuses engendrera des accès mémoires très irréguliers, ce qui dégaradera
540 certainement sa performance de calcul.
541
542 \begin{algorithm}
543   \SetLine
544   \linesnumbered
545   \Entree{$Val$, $Col$ et $Ptr$ (matrice), $n$ (taille de la matrice), $x$ (vecteur)}
546   \Sortie{$y$ (vecteur)}
547   \BlankLine
548    variables $i$, $j$, $id$\; 
549   \Pour{$i=0$ \KwA $n-1$}{
550     $y[i] \leftarrow 0$\;
551   }
552   \BlankLine
553   \Pour{$i=0$ \KwA $n-1$}{
554     \Pour{$id=Ptr[i]$ \KwA $Ptr[i+1]-1$}{
555       $j \leftarrow Col[id]$\;
556       $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
557     }
558   }
559 \caption{Multiplication matrice-vecteur avec le format CSR}
560 \label{alg:csr}
561 \end{algorithm}
562
563 \begin{algorithm}
564   \SetLine
565   \linesnumbered
566   \Entree{$Val$, $Lig$ et $Ptr$ (matrice), $n$ (taille de la matrice), $x$ (vecteur)}
567   \Sortie{$y$ (vecteur)}
568   \BlankLine
569   variables $i$, $j$, $id$\; 
570   \Pour{$i=0$ \KwA $n-1$}{
571     $y[i] \leftarrow 0$\;
572   }
573   \BlankLine
574   \Pour{$j=0$ \KwA $n-1$}{
575     \Pour{$id=Ptr[j]$ \KwA $Ptr[j+1]-1$}{
576       $i \leftarrow Lig[id]$\;
577       $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
578     }
579   }
580 \caption{Multiplication matrice-vecteur avec le format CSC}
581 \label{alg:csc}
582 \end{algorithm}
583
584 \subsection{ELLPACK/ITPACK}
585 \textit{ELLPACK/ITPACK}, ou plus couramment \textit{ELL}, est parmi les structures de données
586 les plus génériques pour le stockage des matrices creuses dans la mémoire. De plus, elle est plus
587 particulièrement adaptée à la représentation des matrices, plus ou moins, structurées et aux
588 calculs sur des architectures vectorielles. Le stockage en format \textit{ELL} d'une matrice creuse,
589 $A$, de taille $(n\times n)$, ayant au maximum $m$ éléments non nuls par ligne, utilise deux matrices
590 de données de taille $(n\times m)$: une matrice de réels, $Val$, et une matrice d'entiers, $Col$.
591 Une ligne d'indice $i$ de chacune des matrices $Val$ et $Col$ correspond aux éléments non nuls et
592 leurs indices de colonne, respectivement, de la ligne d'indice $i$ de la matrice creuse $A$.
593 Les lignes qui ont moins de $m$ éléments non nuls dans la matrice creuse $A$ sont complétées par
594 des zéros dans les matrices $Val$ et $Col$. La figure~\ref{fig:ell} montre un stockage \textit{ELL}
595 d'une matrice creuse de taille $(4\times 4)$. Dans cet exemple, les matrices $Val$ et $Col$ sont
596 complétées par des étoiles ($*$), qui correspondent aux zéros dans la mémoire. L'algorithme~\ref{alg:ell}
597 montre la mutiplication matrice-vecteur
598 avec le format \textit{ELL}, telles que les deux matrices $Val$ et $Col$ soient stockées ligne par ligne
599 dans la mémoire. 
600 \begin{figure}[!h]
601 \begin{center}
602 \begin{tabular}{c}
603 \begin{tabular}{c}
604 $A=\left[
605 \begin{matrix}
606 a_{00} & 0      & 0     & a_{03} \\
607 0     & 0      & a_{12} & 0 \\
608 a_{20} & a_{21} & a_{22} & 0 \\
609 0     & a_{31} & 0      & a_{33} \\
610 \end{matrix}
611 \right]
612 $
613 \end{tabular}
614 \\
615 \\
616 \begin{tabular}{cc}
617 $Val=\left[
618 \begin{matrix}
619 a_{00} & a_{03} & * \\
620 a_{12} & *      & * \\
621 a_{20} & a_{21} & a_{22} \\
622 a_{31} & a_{33} & * \\
623 \end{matrix}
624 \right],
625 $
626 &
627 $Col=\left[
628 \begin{matrix}
629 0 & 3 & * \\
630 2 & * & * \\
631 0 & 1 & 2 \\
632 1 & 3 & * \\
633 \end{matrix}
634 \right]
635 $
636 \end{tabular}
637 \end{tabular}
638 \caption{Le format de stockage ELL}
639 \label{fig:ell}
640 \end{center}
641 \end{figure}
642
643 \begin{algorithm}[H]
644   \SetLine
645   \linesnumbered
646   \Entree{$Val$ et $Col$ (matrice), $n$ (taille de la matrice), $m$ (nombre maximum d'éléments non nuls par ligne), $x$ (vecteur)}
647   \Sortie{$y$ (vecteur)}
648   \BlankLine
649    variables $i$, $j$, $l$, $id$\; 
650   \Pour{$i=0$ \KwA $n-1$}{
651     $y[i] \leftarrow 0$\;
652   }
653   \BlankLine
654   \Pour{$i=0$ \KwA $n-1$}{
655     \Pour{$l=0$ \KwA $m-1$}{
656       $id \leftarrow n \times i + l$\;
657       $j \leftarrow Col[id]$\;
658       $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
659     }
660   }
661 \caption{Multiplication matrice-vecteur avec le format ELL}
662 \label{alg:ell}
663 \end{algorithm}
664
665 \subsection{HYB}
666 Le format hybride (\textit{HYB}) est la combinaison de deux structures de données, qui sont le format 
667 \textit{ELL} et le format \textit{COO}. Bien que le format de stockage \textit{ELL} soint bien adapté aux
668 multiplications matrice-vecteur calculées sur des architectures vectorielles, son efficacité se dégrade
669 rapidement lorsque le nombre d'éléments non nuls varie substantiellemnt d'une ligne de matrice à une
670 autre, ce qui est souvent le cas avec les matrices creuses non structurées. En revanche, la performance
671 des multiplications matrice-vecteur basées sur le format de stockage \textit{COO} est insensible aux
672 variations du nombre d'éléments non nuls par ligne de la matrice creuse.
673
674 Le principe de la structure de données \textit{HYB} consiste à stocker un nombre typique $m$ d'éléments
675 non nuls par ligne en format \textit{ELL} et le reste des éléments non nuls des lignes exeptionnelles
676 en format \textit{COO}. Le nombre typique $m$ est souvent déterminé a priori de telle sorte que la
677 majorité de la matrice creuse ait une structure régulière et ainsi, soit facile à représenter en format
678 \textit{ELL}. De cette manière, le format \textit{HYB} hérite l'efficacité du format \textit{ELL} due
679 à ses accès réguliers à la mémoire et la flexibilité du format \textit{COO} qui est insensible à la
680 structure creuse de la matrice. La figure~\ref{fig:hyb} illustre la structure de données \textit{HYB}
681 d'une matrice creuse $A$ d'ordre $(4\times 4)$, utilisant un format \textit{ELL} avec 2 éléments par
682 ligne ($m=2$). L'algorithme~\ref{alg:hyb} présente le code d'une
683 multiplication matrice-vecteur avec le format de stockage \textit{HYB}.
684
685 \begin{figure}[!h]
686 \begin{center}
687 \begin{tabular}{c}
688
689 \begin{tabular}{c}
690 $A=\left[
691 \begin{matrix}
692 a_{00} & 0      & 0     & a_{03} \\
693 0     & 0      & a_{12} & 0 \\
694 a_{20} & a_{21} & a_{22} & 0 \\
695 0     & a_{31} & 0      & a_{33} \\
696 \end{matrix}
697 \right]
698 $
699 \end{tabular}
700 \\
701 \\
702 \begin{tabular}{cc}
703 $
704 ELL:
705 $
706 &
707 \begin{tabular}{ccc}
708 $
709 Val_{ell}=\left[
710 \begin{matrix}
711 a_{00} & a_{03}\\
712 a_{12} & *    \\
713 a_{20} & a_{21} \\
714 a_{31} & a_{33} \\
715 \end{matrix}
716 \right],
717 $
718 &
719 $
720 Col_{ell}=\left[
721 \begin{matrix}
722 0 & 3\\
723 2 & * \\
724 0 & 1 \\
725 1 & 3 \\
726 \end{matrix}
727 \right]
728 $
729 \end{tabular}
730 \\
731 \\
732 $
733 COO:
734 $
735 &
736 \begin{tabular}{ccc}
737 $Val_{coo}=[a_{22}],$ & $Lig_{coo}=[2],$ & $Col_{coo} = [2]$\\
738 \end{tabular}
739 \\
740 \end{tabular}
741 \end{tabular}
742 \caption{Le format de stockage HYB}
743 \label{fig:hyb}
744 \end{center}
745 \end{figure}
746
747 \begin{algorithm}[H]
748   \SetLine
749   \linesnumbered
750   \Entree{$Val_ {ell}$ et $Col_{ell}$ (matrice ELL), $Val_ {coo}$, $Lig_{coo}$ et $Col_{coo}$ (matrice COO), 
751 $n$ (taille de la matrice ELL), $m$ (nombre maximum d'éléments non nuls par ligne dans la matrice ELL), 
752 $nnz$ (nombre d'éléments non nuls dans la matrice COO), $x$ (vecteur)}
753   \Sortie{$y$ (vecteur)}
754   \BlankLine
755   variables $i$, $j$, $l$, $id$\; 
756   \Pour{$i=0$ \KwA $n-1$}{
757     $y[i] \leftarrow 0$\;
758   }
759   \BlankLine
760   \Pour{$i=0$ \KwA $n-1$}{
761     \Pour{$l=0$ \KwA $m-1$}{
762       $id \leftarrow n \times i + l$\;
763       $j \leftarrow Col_{ell}[id]$\;
764       $y[i] \leftarrow y[i] + Val_{ell}[id] \times x[j]$\;
765     }
766   }
767   \BlankLine
768   \Pour{$id=0$ \KwA $nnz-1$}{
769     $i \leftarrow Lig_{coo}[id]$\;
770     $j \leftarrow Col_{coo}[id]$\;
771     $y[i] \leftarrow y[i] + Val_{coo}[id] \times x[j]$\;
772   }
773 \caption{Multiplication matrice-vecteur avec le format HYB}
774 \label{alg:hyb}
775 \end{algorithm}
776 %%-------------------------------------------------------------------------------------------------------%%
777
778 %%-------------------------------------------------------------------------------------------------------%%
779 \section{Parallélisation des méthodes itératives}
780 \label{sec:parallelisation}
781 L'adaptation des solveurs itératifs aux architectures des calculateurs parallèles est devenue une
782 nécessité impérieuse pour la résolution des systèmes linéaires creux de taille toujours croissante.
783 En effet, les calculateurs parallèles fournissent les ressources et la puissance de calcul nécessaires
784 pour une résolution haute performance de grands systèmes linéaires creux. 
785
786 La parallélisation d'un solveur itératif sur une plateforme parallèle, comportant $p$ processeurs, requiert
787 tout d'abord un partitionnement des données du système linéaire à résoudre sur l'ensemble des $p$ processeurs.
788 Le partitionnement de données consiste à attribuer aux différents processeurs des portions, plus ou moins égales,
789 des vecteurs et des matrices impliqués dans le solveur itératif. Par la suite, tous les processeurs procèdent,
790 en parallèle, à la résolution de leurs parties du système linéaire, en appliquant les itérations d'une même méthode
791 itérative mais sur des données différentes. Toutefois, les calculs locaux des différents processeurs sont liés
792 par des dépendances de données, qui permettent la résolution du système global. A cet effet, des points
793 de synchronisation doivent être mis en place dans un solveur itératif parallèle, lors desquels des données partagées
794 sont échangées entre des processeurs voisins.
795
796 Dans ce document, nous nous intéressons à la parallélisation des méthodes itératives sur des plateformes de calcul
797 parallèle à mémoire distribuée. Donc, lors des points de synchronisation, les dépendances de
798 données sont gérées par des communications à passage de messages entre les processeurs. Les algorithmes
799 parallèles des méthodes itératives peuvent être classifiés selon la nature \textit{synchrone} ou \textit{asynchrone}
800 de leurs itérations ainsi que celle de leurs communications~\cite{ref18}.
801
802 \subsection{Méthodes itératives SISC}
803 Les solveurs itératifs parallèles les plus classiques sont ceux dits à Itérations Synchrones
804 et Communications Synchrones (\textit{SISC}). Le synchronisme des itérations décline du fait que chaque
805 processeur ne peut commencer le calcul de sa nouvelle itération que lorsqu'il reçoit, de la
806 part de tous ses voisins, les données partagées calculées à l'itération précédente. Ainsi, tous
807 les processeurs commencent le calcul de la même itération au même temps et effectuent les échanges
808 de données partagées à la fin de chaque itération par le biais des communications globales synchrones.
809
810 \begin{figure}[!h]
811   \centering
812   \includegraphics[width=170mm,keepaspectratio]{Figures/SISC}
813   \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle SISC avec deux processeurs}
814   \label{fig:SISC}
815 \end{figure}
816
817 La figure~\ref{fig:SISC} montre un schéma d'exécution d'un algorithme parallèle \textit{SISC}.
818 Le synchronisme des itérations des solveurs itératifs parallèles \textit{SISC} impliquent,
819 exactement, le même nombre d'itérations que leurs homologues en séquentiel. Par conséquent,
820 les conditions de leur convergence sont faciles à déterminer, vu qu'elles sont identiques
821 à celles des solveurs itératifs séquentiels. Cependant, l'utilisation des communications
822 synchrones peut pénaliser les performances de ce type d'algorithmes. En fait, les
823 communications synchrones engendrent, souvent, des temps d'inactivité des processeurs
824 (les espaces blancs entre deux itérations calculées par un même processeur dans la
825 figure~\ref{fig:SISC}), dues soit à l'attente pour un processeur (ou des processeurs)
826 qu'il soit prêt pour communiquer ou à la vitesse lente des communications synchrones
827 elles-mêmes dans la plateforme de calcul parallèle.
828
829 \subsection{Méthodes itératives SIAC}
830 Un autre type de solveurs itératifs parallèles a été conçu pour améliorer les performances
831 de résolution sur les plateformes à réseau d'interconnection lent et/ou hétérogène. Ce sont
832 les solveurs à Itérations Synchrones et Communications Asynchrones
833 (\textit{SIAC}). Il permettent de réduire les temps d'inactivité des processeurs entre deux itérations
834 successives, en utilisant des communications asynchrones et cela, tout en maintenant le
835 principe des itérations synchrones. En effet, comme les algorithmes \textit{SISC}, un processeur
836 attend toujours la réception de toutes les données partagées, calculées par ses voisins à
837 l'itération précédente, avant de commencer les calculs de la nouvelle itération. Toutefois,
838 les communications synchrones globales utilisées, pour les échanges de données partagées, dans
839 les algorithmes \textit{SISC} sont remplacées par des envois asynchrones et des réceptions bloquantes.
840
841 \begin{figure}[!h]
842   \centering
843   \includegraphics[width=170mm,keepaspectratio]{Figures/SIAC}
844   \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle SIAC avec deux processeurs}
845   \label{fig:SIAC}
846 \end{figure}
847
848 La figure~\ref{fig:SIAC} montre un exemple de schéma d'exécution d'un solveur itératif parallèle
849 \textit{SIAC}. Les conditions de convergence des solveurs itératifs \textit{SIAC} sont identiques à ceux des solveurs
850 itératifs \textit{SISC} et ainsi, à ceux des solveurs itératifs séquentiels. Cependant, les communications
851 asynchrones permettent d'effectuer des chauvechaments entre les calculs et les échanges de données.
852 En fait, un processeur peut envoyer les données partagées à son voisin dès qu'elles soient 
853 prêtes à être utilisées dans les calculs de la nouvelle itération, ce qui permet de réduire les temps
854 d'attente de réception de données entre deux itérations successives.
855
856 \subsection{Méthodes itératives AIAC}
857 Les méthodes itératives parallèles synchrones de type \textit{SISC} ou \textit{SIAC} ont
858 été, largement, utilisées pour la résolution de grands systèmes linéaires creux. Elles sont, relativement,
859 faciles à mettre en \oe uvre et présentent des taux de convergence, plus ou moins, rapides selon la méthode utilisée. Néanmoins, le synchronisme
860 de leurs itérations et/ou de leurs communications pénalise sévèrement les performances de
861 résolution sur des grilles de calcul à clusters géographiquement distants. En effet, ce
862 type d'architectures parallèles est souvent caractérisé par la latence forte de ses liens de
863 communications (plus précisément les liens inter-clusters) et par l'hétérogénéité de ses
864 ressources. 
865
866 Les solveurs itératifs parallèles conçus pour améliorer les performances de résolution
867 sur les plateformes à ressources hétérogènes et géographiquement distantes sont ceux
868 dits à Itérations Asynchrones et Communications Asynchrones (\textit{AIAC}). Dans ce
869 type de solveurs parallèles, chaque processeur exécute ses propres itérations sans
870 prendre en compte la progression de l'exécution de celles des autres processeurs. En
871 effet, un processeur passe d'une itération à l'autre sans attendre l'arrivée des nouvelles
872 données partagées mises à jour par ses voisins. Il utilise les données locales et les
873 dernières versions des données partagées disponibles au début du calcul de chaque itération.
874 De plus, pour assurer l'asynchronisme des itérations, les solveurs itératifs \textit{AIAC}
875 utilisent des communications asynchrones de type envois et réceptions non bloquants.
876
877 \begin{figure}[!h]
878   \centering
879   \includegraphics[width=170mm,keepaspectratio]{Figures/AIAC}
880   \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle AIAC avec deux processeurs}
881   \label{fig:AIAC}
882 \end{figure}
883
884 La figure~\ref{fig:AIAC} présente un schéma d'exécution d'un solveur itératif
885 parallèle de type \textit{AIAC}. Nous pouvons remarquer que l'asynchronisme des
886 itérations et des communications permet aux processeurs d'éviter les temps
887 d'inactivité et d'exécuter des itérations différentes à un moment donné de la
888 résolution. De ce fait, certains processeurs peuvent être plus rapides dans leurs
889 calculs et effectuer plus d'itérations que les autres. Cependant, les solveurs
890 itératifs asynchrones \textit{AIAC} ont des conditions de convergence plus strictes
891 et effectuent plus d'itérations que les solveurs synchrones \textit{SISC} et \textit{SIAC}.
892 Par conséquent, ils nécessitent une analyse plus élaborée pour déterminer les
893 bons indicateurs de convergence.
894 %%-------------------------------------------------------------------------------------------------------%%
895
896 %%-------------------------------------------------------------------------------------------------------%%
897 \section{Conclusion}
898 Dans ce chapitre, nous avons commencé par présenter les méthodes de résolution des systèmes linéaires creux.
899 Premièrement, nous avons abordé l'analyse des méthodes directes qui permettent de trouver la solution exacte en un
900 nombre d'opérations élémentaires fini. Cependant, pour la résolution des systèmes linéaires creux, nous avons mis
901 en évidence l'existence d'une opération de remplissage à l'étape de factorisation des matrices creuses. L'opération de
902 remplissage est très coûteuse en termes de temps d'exécution et d'espace mémoire et, plus précisément, pour
903 les matrices creuses de grande taille. Ensuite, nous avons décrit les méthodes itératives qui permettent de 
904 calculer une solution approximative en exécutant des itérations successives d'un même bloc d'opérations
905 élémentaires. Nous avons cité trois grandes familles de méthodes itératives à savoir: les méthodes 
906 stationnaires, les méthodes non stationnaires et les méthodes multigrilles. Celles-ci sont plus adaptées à la
907 résolution des systèmes linéaires de grande taille et plus faciles à paralléliser. 
908
909 Puis, nous avons étudié différentes structures de données des matrices creuses en mémoire. Elles permettent
910 un stockage optimisé des éléments non nuls des matrices creuses, de façon à ce que les accès en lecture/écriture
911 à la mémoire soient faciles à effectuer. Par conséquence, elles permettent d'améliorer les performances de
912 calcul des méthodes de résolution. Enfin, nous avons donné les principaux points clés de la parallélisation
913 des méthodes itératives sur les plateformes de calcul parallèle à mémoire distribuée, afin de résoudre des
914 systèmes linéaires creux à très grande échelle. Nous avons présenté trois différentes méthodes itératives
915 parallèles, classifiées selon la nature synchrone ou asynchrone de leurs itérations et/ou de leurs communications.
916 Les méthodes synchrones sont plus efficaces sur des petits clusters homogènes ayant des liens de communications
917 rapides, alors que les méthodes asynchrones offrent de meilleures performances sur les plateformes à clusters
918 hétérogènes et géographiquement distants.
919
920 Dans la suite de ce document, nous nous intéresserons à la parallélisation des méthodes itératives sur des
921 clusters équipés de cartes graphiques GPUs. Nous décrirons les principaux points clés des algorithmes parallèles
922 de quelques méthodes itératives, les plus utilisées, pour la résolution des systèmes linéaires ou non linéaires
923 creux. Ensuite, nous effectuerons une comparaison de performance entre les solveurs itératifs parallèles
924 mis en \oe uvre sur les clusters de GPUs et les mêmes solveurs mis en \oe uvre sur les clusters CPUs.      
925 %%-------------------------------------------------------------------------------------------------------%%
926
927 %%% Local Variables: 
928 %%% mode: latex
929 %%% TeX-master: "these"
930 %%% End: