\chapter{R\'esolution de syst\`emes lin\'eaires creux}
\label{chap:sys_lin_creux}

%%-------------------------------------------------------------------------------------------------------%%
%Introduction
Les systèmes linéaires creux sont utilisés pour la modélisation de nombreux problèmes scientifiques ou industriels,
comme la simulation des phénomènes de l'environnement et de la terre ou la transformation industrielle des
fluides complexes ou non-newtoniens. De plus, la résolution de ces problèmes fait souvent appel à la résolution
de ces systèmes linéaires, qui est considérée comme le processus le plus coûteux en termes de temps de calcul
et d'espace mémoire. Par conséquent, la résolution de systèmes linéaires creux doit être aussi efficace que
possible, afin de pouvoir traiter des problèmes de taille toujours croissante.

Il existe, dans le jargon de l'analyse numérique, différentes méthodes de résolution des systèmes linéaires
creux. Cependant, le choix d'une méthode de résolution est souvent guidé par les propriétés de la matrice creuse
du système à résoudre (structure, symétrie, densité, dimension, etc), la vitesse et la précision
de résolution souhaitées. Toutefois, un bon choix de la méthode est aussi guidé par la taille de ces systèmes,
qui ralentit davantage le processus de résolution. Par conséquence, ces méthodes doivent être faciles à paralléliser
pour exploiter la puissance de calcul et la capacité mémoire des plateformes de calcul parallèle.  

Dans le présent chapitre, nous décrivons deux classes de méthodes de résolution dans
la section~\ref{sec:methodes_resolu}, à savoir les méthodes directes et les méthodes itératives. Dans la
section~\ref{sec:stockage}, nous montrons les différentes structures de données des matrices creuses
en mémoire, qui permettent d'optimiser l'occupation d'espace mémoire et les accès en lecture/écriture
à la mémoire. Enfin dans la section~\ref{sec:parallelisation}, nous présentons la parallélisation des 
méthodes itératives pour la résolution de systèmes linéaires creux de grande taille sur les plateformes de 
calcul parallèle.
%%-------------------------------------------------------------------------------------------------------%%

%%-------------------------------------------------------------------------------------------------------%%
\section{Méthodes de résolution}
\label{sec:methodes_resolu}
Soit le système linéaire creux suivant:
\begin{equation}
Ax=b,
\label{eq:linear_syst}
\end{equation}
où $A\in\mathbb{R}^{n\times n}$ est une matrice carrée creuse inversible, $x\in\mathbb{R}^{n}$ est le vecteur
solution, $b\in\mathbb{R}^{n}$ est le vecteur du second membre et $n\in\mathbb{N}$ est un grand entier
naturel. Les méthodes de résolution des systèmes linéaires sont, plus ou moins, classifiées en deux grandes
familles: les méthodes \emph{directes} et les méthodes \emph{itératives}. Dans l'analyse numérique,
un algorithme permettant de résoudre un système linéaire creux est appelé un \emph{solveur linéaire creux direct}
ou \emph{itératif}, selon le type de la méthode de résolution qu'il utilise.

\subsection{Méthodes directes}
Les méthodes directes sont les méthodes numériques les plus anciennes pour la résolution des systèmes
linéaires. De plus, jusqu'aux années quatre-vingts,
%dans les années précédentes, 
elles ont souvent été préférées aux méthodes itératives,
en raison de leur \emph{robustesse} et de leur comportement \emph{prévisible}. En effet, que ce soit
pour des matrices denses ou creuses, ces méthodes permettent d'obtenir, théoriquement en l'absence des erreurs
d'arrondi, une solution \emph{exacte} en un nombre d'opérations élémentaires \textit{fini} et \emph{connu a priori}.
La plupart des méthodes directes sont basées, principalement, sur deux étapes de calcul: la factorisation
de la matrice suivant ses propriétés (symétrie, définie positive, etc) et la résolution
par des substitutions successives~\cite{ref1}. 

Parmi les solveurs directs les plus utilisés, nous pouvons citer le \textit{solveur LU} basé sur la
méthode de factorisation \textit{LU}~\cite{ref2}. Dans un premier temps, ce solveur procède à la factorisation de
la matrice des coefficients $A$ du système linéaire à résoudre (\ref{eq:linear_syst}). Dans cette étape,
la matrice $A$ est décomposée en deux matrices triangulaires inférieure $L=(l_{ij})_{1\leqslant i,j\leqslant n}$
et supérieure $U=(u_{ij})_{1\leqslant i,j\leqslant n}$, de façon à ce que la matrice $A$ puisse s'écrire sous
forme d'un produit des deux matrices $L$ et $U$. Une fois la factorisation établie, le solveur \textit{LU}
résout le système linéaire en deux étapes:
\begin{itemize*}
\item la résolution du système $Ly=b$ par l'algorithme de \emph{descente} (substitutions successives
avant ou résolution de haut vers le bas),
\item puis, la résolution du système $Ux=y$ par l'algorithme de \emph{remontée} (substitutions successives
arrière ou résolution de bas vers le haut).
\end{itemize*}

Cependant, les méthodes directes de résolution des systèmes linéaires creux sont souvent plus
compliquées que celles des systèmes denses. La principale difficulté réside dans l'étape de factorisation
de la matrice. En effet, durant la construction des deux facteurs triangulaires $L$ et $U$, certains
zéros de la matrice creuse d'origine, $A$, peuvent devenir des valeurs non nulles dans les matrices $L$
et $U$, ce que l'on appelle, communément, processus de \emph{remplissage}. Par conséquent, une mise en \oe uvre
efficace d'un solveur direct creux dépend largement des structures de données en mémoire (formats de stockage)
des valeurs non nulles des matrices $L$ et $U$. De ce fait, un solveur direct creux est généralement
basé sur quatre étapes distinctes~\cite{ref3,ref4} au lieu de deux comme dans les solveurs denses:
\begin{enumerate*}
 \item Réorganisation de la matrice $A$: qui consiste à réordonner les lignes et les colonnes de $A$,
et donc les valeurs non nulles, de telle sorte que les facteurs $L$ et $U$ subissent peu de remplissage
à l'étape de factorisation numérique;
 \item Factorisation symbolique: qui détermine a priori les structures de données appropriées aux valeurs
non nulles des facteurs $L$ et $U$;
 \item Factorisation numérique: qui calcule les facteurs $L$ et $U$;
 \item Résolution du système linéaire creux: qui effectue les substitutions avant-arrière en utilisant
les deux facteurs $L$ et $U$ calculés à l'étape précédente. 
\end{enumerate*}
Dans le cas des matrices creuses symétriques et définies positives (valeurs propres strictement posistives),
les quatres étapes peuvent être bien
séparées comme dans le solveur \textit{CHOLMOD}~\cite{ref5}. Pour les matrices creuses asymétriques,
les solveurs directs peuvent fusionner les étapes 2 et 3 en une seule étape comme les solveurs \textit{SuperLU}~\cite{ref6}
et \textit{MUMPS}~\cite{ref20},
ou même fusionner les étapes 1, 2 et 3 comme le solveur \textit{UMFPACK}~\cite{ref7}.

Toutefois, les méthodes directes ne sont pas bien adaptées pour la résolution des systèmes linéaires
creux de très grandes dimensions, et en particulier, ceux associés aux problèmes tridimensionnels. En fait,
l'occupation espace mémoire du processus de factorisation croît \textit{superlinéairement}
avec la taille de la matrice et dépend fortement de la structure de cette dernière (positions des valeurs non nulles)
pour réduire le remplissage. Hélas, quelque soit le nombre de fois que les coefficients de la
matrice sont réordonnés, cela n'empêche pas le remplissage de croître rapidement et proportionnellement
à la taille de la matrice. En outre, les solveurs directs sont très coûteux en termes de temps de résolution
qui est de $\mathcal{O}(n^{3})$, où $n$ représente la taille du système. Dans ce cas, les méthodes itératives
décrites dans la section~\ref{sec:iterative}, sont particulièrement plus utiles lorsque les systèmes linéaires
à résoudre sont creux et de très grande dimension car, contrairement aux méthodes directes, aucun processus de remplissage n'est nécessaire.  

\subsection{Méthodes itératives}
\label{sec:iterative}
Les méthodes itératives procédent par itérations successives
d'un même bloc d'opérations élémentaires, au cours desquelles une séquence infinie de solutions approximatives
$\{x^{(k)}\}_{k\geq0}$ est générée. A partir d'une solution initiale $x^{(0)}$ (choisie pour la première
itération), une méthode itérative détermine à chaque itération
$k\geq1$, une solution approximative, $x^{(k)}$, qui converge progressivement vers la solution exacte, $x^{*}$,
du système linéaire~\ref{eq:linear_syst}, telle que:
\begin{equation}
x^{*}=\lim\limits_{k\to\infty}x^{(k)}=A^{-1}b\in\mathbb{R}^{n}.
\end{equation}

A la différence d'une méthode directe, le nombre d'itérations d'une méthode itérative nécessaire pour
atteindre la solution exacte du système linéaire à résoudre \textit{n'est pas connu} au préalable et
peut être \textit{infini}. Cependant, dans la pratique, un solveur itératif donne souvent une solution
approximative satisfaisant une précision requise et/ou après un nombre maximum d'itérations fixé a priori.
Les critères d'arrêt les plus fréquemment utilisés par les solveurs itératifs sont les suivants:
\begin{itemize}
  \item $\|x^{(k)}-x^{(k-1)}\|<\varepsilon$: erreur absolue,
  \item $\|x^{(k)}-x^{(k-1)}\|<\varepsilon\|x^{(k-1)}\|$: erreur relative,
  \item $\|b-Ax^{(k)}\|<\varepsilon\|b\|$: résidu de la solution. 
\end{itemize}
où $x^{(k)}$ est la solution approximative obtenue à l'itération $k\leqslant K$, $K$ étant le nombre maximum
d'itérations, et $\varepsilon$ est le seuil de tolérance requis pour la précision de calcul. Ainsi, un solveur
linéaire itératif est dit \textit{convergent} si pour tout choix de la solution initiale $x^{(0)}\in\mathbb{R}^{n}$,
l'erreur ou le résidu de la solution finale $x^{(k)}\in\mathbb{R}^{n}$, obtenue après $k$ itérations, est
suffisamment \textit{petit}. Cependant, la convergence des méthodes itératives n'est pas toujours assurée
et dépend fortement des caractéristiques du système linéaire, à savoir: la structure creuse, le conditionnement
et le rayon spectral de la matrice des coefficients $A$. Il existe différentes classes de méthodes itératives
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
linéaire à résoudre. Nous citons ici, globalement, trois classes de méthodes itératives:
stationnaires, non stationnaires et multigrilles.

\subsubsection{Méthodes stationnaires}
Les méthodes itératives \textit{stationnaires} permettant de résoudre le système linéaire (\ref{eq:linear_syst}),
peuvent être exprimées sous la forme simple suivante:
\begin{equation}
\forall k\in\mathbb{N}, x^{(k+1)}=f(x^{(k)})=Cx^{(k)}+d,
\end{equation}
telle que $f$ soit une fonction \textit{affine} représentant un schéma de point fixe, où la matrice d'itération $C$
et le vecteur $d$ restent \textit{inchangés} pour toutes les itérations $k$. Les méthodes itératives
stationnaires sont basées sur la décomposition de la matrice des coefficients, $A$, comme suit:
\begin{equation}
A=M-N,
\end{equation}
de sorte que le système linéaire (\ref{eq:linear_syst}) puisse s'écrire sous forme: $Mx=Nx+b$ à laquelle est 
associée l'itération suivante:
\begin{equation}
x^{(k+1)}=Cx^{(k)}+d=M^{-1}Nx^{(k)}+M^{-1}b,
\end{equation}
où $M$ est une matrice inversible. Par conséquent, une méthode itérative stationnaire est \textit{convergente}
vers la solution exacte $x^{*}$ si, et seulement si, le rayon spectral (la plus grande valeur propre) de la matrice $C$ est strictement
inférieur à $1$: $\rho(C)<1$. Les principales méthodes itératives stationnaires sont basées sur la décomposition
de la matrice des coefficients en matrices diagonale $D$, triangulaire strictement inférieure $L$ et
triangulaire strictement supérieure $U$: $A=D-L-U$. Parmi ces méthodes itératives nous pouvons citer~\cite{ref8}:
\begin{itemize}
  \item Jacobi: $M=D$ et $N=L+U$,
  \item Gauss-Seidel: $M=D-L$ et $N=U$, 
  \item Surrelaxation successive (SOR): $M=D-\omega L$ et $N=\omega U+(1-\omega)D$ et $\omega>0$, 
  \item Richardson: $C=I-\omega A$, $d=\omega b$ et $\omega>0$.
\end{itemize}

L'algorithme~\ref{alg:stationnaire} présente les principaux points clés d'un solveur itératif
stationnaire. Il permet de calculer, à partir d'une solution initiale $x^{(0)}$, une série de
solutions approximatives successives $x$, jusqu'à ce que la convergence soit atteinte
(dans l'algorithme, à la ligne $7$, résidu<$\varepsilon$).

\begin{algorithm}[H]
  \SetLine
  \linesnumbered
  \Entree{$A$ (matrice), $C$ (matrice d'itération), $b$ (vecteur), $d$ (vecteur)}
  \Sortie{$x$ (vecteur)}
     
  \BlankLine

  choisir $x^{(0)}$ solution initiale, $\varepsilon$ seuil de tolérance\;
  variables $k$, $convergence$\;
  $convergence \leftarrow faux$\;
  $k \leftarrow 0$\;
  
  \Tq{$\neg convergence$}
  {
    évaluer: $x^{(k+1)}=f(x^{(k)})=Cx^{(k)}+d$\;
    $convergence \leftarrow (\|b-Ax^{(k+1)}\|\cdot\|b\|^{-1}<\varepsilon)$\;
    $k \leftarrow k+1$\;
  }
\caption{Algorithme général d'un solveur itératif stationnaire}
\label{alg:stationnaire}
\end{algorithm}

La ligne $6$ de l'algorithme~\ref{alg:stationnaire} permet de calculer les composantes du
vecteur $x$ en appliquant la fonction affine de la méthode itérative stationnaire utilisée.
Par exemple, dans les méthodes Jacobi et Gauss-Seidel, les $n$ composantes du vecteur
$x$ à l'itération $(k+1)$ sont calculées comme suit: 

\begin{tabular}{rl}
  $Jacobi:$& $ x^{(k+1)}=D^{-1}(L+U)x^{(k)}+D^{-1}b$\hspace{1.6cm}$\Leftrightarrow$\\ \\
  & $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$,  \\ \\
  $Gauss-Seidel:$ & $x^{(k+1)}=(D-L)^{-1}Ux^{(k)}+(D-L)^{-1}b$\hspace{0.5cm}$\Leftrightarrow$\\ \\
  & $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$. 
\end{tabular}

Comme nous l'avons évoqué, précédemment, le principal critère de choix entre les solveurs
itératifs est l'efficacité, et donc la convergence, des méthodes de résolution utilisées
par ces derniers. Par exemple, pour trouver la nouvelle solution $x$ à l'itération $(k+1)$,
la méthode Jacobi utilise toutes les composantes du vecteur $x$ calculées à l'itération
précédente $k$, alors que la méthode Gauss-Seidel utilise les nouvelles composantes
calculées à l'itération courante $(k+1)$ dès que possible. Par conséquent, cette technique
de mise à jour des composantes du vecteur solution $x$ permet au solveur Gauss-Seidel
d'effectuer moins d'itérations et donc d'avoir, en général, une convergence plus rapide que
le solveur Jacobi.

Les méthodes stationnaires sont les méthodes itératives les plus anciennes et leurs
algorithmes sont simples à comprendre et à mettre en \oe uvre. Cependant, elles sont souvent
moins efficaces pour la résolution de beaucoup de systèmes linéaires creux, vu leur applicabilité
limitée à certains types de matrices. Ces dernières années, les méthodes non stationnaires,
définies dans la section suivante, sont devenues les méthodes itératives les plus utilisées
pour leur efficacité à résoudre plusieurs types de systèmes linéaires.

\subsubsection{Méthodes non stationnaires}
\label{sec:non-stationnaire}
Les méthodes \textit{non stationnaires} sont les méthodes itératives les 
plus récentes; leurs algorithmes sont généralement plus difficiles à comprendre et à 
mettre en \oe uvre, mais elles peuvent être plus efficaces pour la résolution de certains systèmes
linéaires creux. Elles peuvent être exprimées sous la forme suivante:
\begin{equation}
  \forall k\in\mathbb{N}, x^{(k+1)}=f_{k}(x^{(k)}),
\end{equation}
où, contrairement aux méthodes stationnaires, la fonction $f_{k}$ permettant de calculer
le nouvel itéré $x^{(k+1)}$ diffère d'une itération à l'autre. En fait, les méthodes non stationnaires 
impliquent des données qui \textit{changent} à chaque itération et qui déclinent des opérations
vectorielles utilisées par les méthodes itératives, comme les produits scalaires et les 
opérations \textit{AXPY} ($y\leftarrow ax+y$).

Le principe général des méthodes non stationnaires est la construction d'une séquence de 
vecteurs orthogonaux et les projections sur des sous-espaces. L'une des classes importantes
des méthodes non stationnaires est celle des \textit{méthodes itératives de Krylov}.
L'idée générale d'une méthode de Krylov~\cite{ref8} est de chercher une approximation, $x^{(k)}$, de 
la solution de système linéaire (\ref{eq:linear_syst}) dans un sous-espace de Krylov,
$\mathcal{K}_{k}$, d'ordre $k$:
\begin{equation}
x^{(k)}-x^{(0)}=q_{k-1}(A)r^{(0)}\in\mathcal{K}_{k}(A,r^{(0)}),
\end{equation}
telle que la condition Petrov-Galerkin soit satisfaite:\[b-Ax^{(k)}\bot\mathcal{L}_{k},\]
où $q_{k-1}$ est un polynôme de degré $k-1$ et $\mathcal{L}_{k}$ est un autre sous-espace
vectoriel d'ordre $k$. Le sous-espace de Krylov d'ordre $k$, noté $\mathcal{K}_{k}(A,r^{(0)})$, 
est un sous-espace vectoriel construit sur $A$ et $r^{(0)}$ et engendré par les $k$ vecteurs
orthogonaux deux à deux $r^{(0)}$, $Ar^{(0)}$, $\dots$, $A^{(k-1)}r^{(0)}$ (base orthonormée):
\[\mathcal{K}_{k}(A,r^{(0)})\equiv \text{engendré}\{r^{(0)}, Ar^{(0)}, A^{2}r^{(0)}, \dots, A^{k-1}r^{(0)}\},\]
où $r^{(0)}=b-Ax^{(0)}$ est le résidu initial de la solution initiale $x^{(0)}$ choisie, plus 
ou moins aléatoirement, au démarrage du solveur de Krylov. 

La plupart des méthodes de Krylov sont basées sur des itérations ne requierant que des multiplications
de la matrice du système par un vecteur, des produits scalaires et des additions de vecteurs; et
elles diffèrent généralement dans le choix de $\mathcal{L}_{k}$, le calcul d'une base orthonormée
pour $\mathcal{K}_{k}$ (certaines méthodes utilisent la matrice creuse $A$ et d'autres utilisent 
sa transposée $A^{T}$) et la réinitialisation de cette base à toutes les $m$ itérations (restarts).
Les méthodes de Krylov ont, généralement, une convergence plus rapide que les méthodes stationnaires.
Néanmoins, pour certains types de matrices creuses, elles n'arrivent pas à converger ou elles
effectuent de nombreuses itérations pour converger. En pratique, elles sont, généralement,
utilisées avec un \textit{préconditionneur} qui permet d'améliorer et/ou d'accélérer leur convergence.
Le processus de préconditionnement consiste à remplacer le système linéaire creux (\ref{eq:linear_syst})
par l'un des deux systèmes suivants:
\begin{equation}M^{-1}Ax=M^{-1}b,\label{eq:left}\end{equation}
\begin{center}ou\\\end{center}
\begin{equation}\begin{array}{cc}AM^{-1}\hat x=b, & x=M^{-1}\hat x,\end{array}\label{eq:right}\end{equation}
qui sont plus faciles à résoudre. $M$ représente la matrice de préconditionnement et les formules (\ref{eq:left})
et (\ref{eq:right}) représentent, respectivement, le préconditionnement à gauche et le préconditionnement
à droite du système linéaire $Ax=b$. Il existe beaucoup de méthodes de Krylov~\cite{ref8,ref9}, par exemple
\textit{CG} ($\mathcal{L}_{k}=\mathcal{K}_{k}$), \textit{GMRES} ($\mathcal{L}_{k}=A\mathcal{K}_{k}$) et \textit{BICG}
($\mathcal{L}_{k}=\mathcal{K}_{k}(A^{T},r^{(0)})$).

\subsubsection{Méthodes multigrilles}
Une autre catégorie de méthodes itératives est celle des méthodes \textit{multigrilles} dont les origines
remontent aux années soixante \cite{ref10}. Initialement conçues pour la résolution de l'équation de Poisson,
les méthodes multigrilles ont rapidement évolué pour être appliquées à un large éventail d'applications dans
le domaine de la physique ou des mathématiques des finances, comme par exemple la résolution des problèmes
linéaires ou non linéaires avec des conditions aux limites \cite{ref11} ou des problèmes d'optimisation \cite{ref12}.

Les méthodes itératives standards, comme Jacobi et Gauss-Seidel, sont caractérisées par un faible taux de
convergence global, $\tau$, qui peut être défini à l'itération $k$ comme suit:
\begin{equation}
  \tau_{k}=\frac {\|r^{(k)}\|} {\|r^{(k-1)}\|} = \frac {\|b-Ax^{(k)}\|} {\|b-Ax^{(k-1)}\|},
\end{equation}
où $r^{(k)}$ est le résidu de la solution approximative $x^{(k)}$ trouvée après $k$ itérations. En fait, malgré
leur convergence rapide aux premières itérations, ces méthodes itératives peinent à maintenir leur vitesse
de convergence tout le long de la résolution. Cela peut être expliqué par leur incapacité à réduire les basses
fréquences de l'erreur, ce qui se traduit par la décroissance lente de la norme du résidu durant les dernières
itérations de la résolution. Par conséquent, beaucoup d'itérations sont nécessaires pour donner une solution
approchée intéressante lorsque la vitesse de convergence est faible. A cet effet, les méthodes multigrilles
sont conçues pour améliorer la vitesse de convergence de ces méthodes itératives et ainsi, réduire leur coût
de résolution. 

Comme leur nom l'indique, les méthodes multigrilles utilisent différents niveaux de maillages réguliers
(au moins deux grilles), du plus fin au plus grossier, chaque n\oe ud de maillage étant espacé d'un autre
d'une distance $2qh$, tel que $q\in\mathbb{N}$ représente le niveaux de la grille. Pour des raisons de
simplicité, nous ne décrirons, ci-après, que le principe général de résolution d'une méthode à deux grilles
sur laquelle sont basées, d'ailleurs, l'ensemble des méthodes multigrilles. On considère un système
linéaire creux (\ref{eq:fine}) discrétisé sur un maillage $\Omega_{h}$ de pas $h=\frac{1}{n+1}$, $n$
étant le nombre de points de discrétisation:
\begin{equation}
  \begin{array}{lll}
    A_{h}x_{h}=b_{h}, & sur & \Omega_{h},\\
  \end{array} 
\label{eq:fine}
\end{equation}
où $A_{h}$, $x_{h}$ et $b_{h}$ représentent, respectivement, la matrice creuse, le vecteur solution et le
second membre du système discrétisé sur $\Omega_{h}$. La résolution de ce système par une méthode itérative
permet de trouver seulement une approximation $x_{h}$ à la solution exacte $x_{h}^{*}$ sur $\Omega_{h}$.
L'erreur algébrique $e_{h}$ et le résidu $r_{h}$ de cette solution approximative peuvent être calculés comme
suit: 
\begin{equation}
  e_{h}=x^{*}_{h}-x_{h},
\label{eq:erreur}
\end{equation}
\begin{equation}
  r_{h}=b_{h}-A_{h}x_{h},
\label{eq:residu}
\end{equation}
De ces deux dernières formules, nous pouvons déduire l'équation de l'erreur suivante:
\begin{equation}
 \begin{array}{lll}
  A_{h}e_{h}=r_{h}, & sur & \Omega_{h},
 \end{array} 
\label{eq:residuelle}
\end{equation} 
L'objectif de la méthode bigrille \cite{ref13} est de réduire les basses fréquences de l'erreur sur
la grille grossière $\Omega_{2h}$, que la méthode itérative n'a pas réussie à éliminer sur la gille fine
$\Omega_{h}$. En fait, les basses fréquences de l'erreur sur la grille fine apparaissent comme des hautes fréquences de
l'erreur sur la grille grossière, qui sont faciles à éliminer par les méthodes itératives. Le principe
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
grossière $\Omega_{2h}$, pour trouver une erreur $e_{h}$ permettant de corriger la valeur de la solution
approximative $x_{h}$: $x_{h}=x_{h}+e_{h}$, obtenue après quelques itérations sur la grille fine
$\Omega_{h}$. L'algorithme \ref{alg:bigrille} montre les principaux points clés d'une méthode bigrille.
Les processus de pré-lissage et de post-lissage consistent à effectuer quelques itérations d'un solveur
itératif (par exemple Jacobi ou Gauss-Seidel) sur la grille fine, afin d'éliminer les hautes fréquences
de l'erreur. Les deux opérateurs de restriction $R^{2h}_{h}$ et de prolongement $P^{h}_{2h}$ permettent le
transfert des données entre la grille fine et la grille grossière. $A_{2h}$ est la matrice d'itération sur
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
avec une méthode itérative sur la grille grossière (à la ligne $7$ de l'algorithme \ref{alg:bigrille})
permet de trouver l'erreur de correction de la solution approximative $x_{h}$ obtenue sur la grille fine
et ainsi, d'éliminer les basses fréquences de l'erreur. 

\begin{algorithm}[H]
  \SetLine
  \linesnumbered
  \Entree{$A$ (matrice), $b$ (vecteur)}
  \Sortie{$x$ (vecteur)}
     
  \BlankLine

  choisir $x_{h}^{0}$ solution initiale, $\varepsilon$ seuil de tolérance\;
  $convergence \leftarrow faux$\;
  
  \Tq{$\neg convergence$}
  {
    Pré-lissage: résoudre $A_{h}x_{h}=b_{h}$ sur $\Omega_{h}$\;
    Calcul du résidu: $r_{h}=b_{h}-A_{h}x_{h}$\;
    Restriction sur $\Omega_{2h}$: $r_{2h}=R_{h}^{2h}r_{h}$\;
    Résoudre: $A_{2h}e_{2h}=r_{2h}$ sur $\Omega_{2h}$\;
    Prolongement sur $\Omega_{h}$: $e_{h}=P^{h}_{2h}e_{2h}$\;
    Correction: $x_{h}=x_{h}+e_{h}$\;
    Post-lissage: résoudre $A_{h}x_{h}=b_{h}$ sur $\Omega_{h}$\; 
    $convergence \leftarrow (\|b_{h}-A_{h}x_{h}\|<\varepsilon)$\;
    $x_{h}^{0}\leftarrow x_{h}$\;
  }
\caption{Algorithme général d'une méthode bigrille}
\label{alg:bigrille}
\end{algorithm}

Cependant, dans la pratique, les méthodes multigrilles ne se limitent pas à deux grilles
pour obtenir une convergence rapide. Elles sont, souvent, appliquées avec plusieurs grilles
de plus en plus grossières: $\Omega_{h}$, $\Omega_{2h}$, $\Omega_{4h}$, $\Omega_{8h}$, $\cdots$.
L'idée générale de ces méthodes est l'application récursive du principe de la méthode bigrille.
En effet, les méthodes multigrilles appliquent une méthode bigrille pour résoudre l'équation
de l'erreur (ligne 7 de l'algorithme~\ref{alg:bigrille}), afin de trouver une meilleure correction
de la solution approximative. Il existe plusieurs versions de ce type de méthodes, par exemple
\textit{V-cycle}, \textit{W-cycle} ou \textit{FMG}, qui sont présentées dans~\cite{ref14,ref15,ref16}.
%%-------------------------------------------------------------------------------------------------------%%

%%-------------------------------------------------------------------------------------------------------%%
\section{Formats de stockage des matrices creuses}
\label{sec:stockage}
La plupart des méthodes itératives sont basées sur des opérations de calcul traitant des vecteurs et des
matrices (voir section~\ref{sec:iterative}). Cependant, parmi toutes ces opérations, seule la multiplication
matrice-vecteur est considérée comme l'opération la plus importante des méthodes itératives, car elle est
souvent très coûteuse en termes de temps d'exécution et d'occupation espace mémoire. Par conséquent,
l'efficacité d'une méthode itérative dépend fortement de la performance de calcul de sa multiplication
matrice-vecteur.

En outre, dans le cas des systèmes linéaires creux, la multiplication matrice-vecteur nécessite, aussi, une
attention très particulière pour le format de stockage de la matrice creuse dans la mémoire. Le stockage
naïf, ligne par ligne ou colonne par colonne, d'une matrice creuse peut causer une perte considérable
d'espace mémoire et de temps de calcul. En fait, les éléments nuls de la matrice sont sauvegardés dans
la mémoire alors que leur prise en compte dans la multiplication est totalement inutile. De plus, la
nature creuse de la matrice engendre, souvent, des accès irréguliers à la mémoire pour la lecture de ses
éléments non nuls et ceci, ne peut que dégrader encore davantage la performance de la multiplication
matrice-vecteur.

Néanmoins, une matrice creuse permet, aussi, de tirer profit de son aspect creux et ce, en ne sauvegardant
que ses éléments non nuls, qui permettra ainsi, de réduire la taille de la mémoire nécessaire
pour le stockage et d'accélérer le calcul de la multiplication matrice-vecteur. Il existe dans la
littérature~\cite{ref8,ref17}, plusieurs formats de stockage compressés pour les matrices creuses. Ils permettent un accès
facile aux éléments non nuls en connaissant leurs coordonnées dans la mémoire et une exécution intelligente
des opérations de calcul sur ces éléments. Dans la présente section, nous décrivons au total cinq
formats de compression des matrices creuses, qui sont intéressants pour la suite de ce document.

\subsection{COO}
\textit{COO} (Coordinate) est un format de stockage particulièrement simple. Il permet de représenter une matrice
creuse sous forme d'une structure de données composée de trois tableaux: un tableau de réels, $Val$, contenant
tous les éléments non nuls de la matrice creuse $A$ et deux tableaux d'entiers,
$Lig$ et $Col$, contenant les indices de lignes et les indices de colonnes de ces éléments non nuls, respectivement.
Les trois tableaux ont une taille $nnz$ (nombre d'éléments non nuls de la matrice). La figure~\ref{fig:coo}
décrit le format \textit{COO} pour une matrice creuse d'ordre $(4\times4)$, $A$, et l'algorithme~\ref{alg:coo} présente
la multiplication matrice-vecteur pour ce format de stockage. 
\begin{figure}[!h]
\begin{center}
\begin{tabular}{c}
$A=\left[
\begin{matrix}
a_{00} & 0 & 0 & a_{03} \\
0 & 0 & a_{12} & 0 \\
a_{20} & a_{21} & a_{22} & 0 \\
0 & a_{31} & 0 & a_{33} \\
\end{matrix}
\right]
$
\\
\\
$
\begin{array}{cccccccccccc}
Val & = & [ & a_{00} & a_{03} & a_{12} & a_{20} & a_{21} & a_{22} & a_{31} & a_{33} & ] \\
Lig & = & [ & 0     & 0      & 1     & 2      & 2     & 2      & 3     & 3      & ] \\
Col & = & [ & 0     & 3      & 2     & 0      & 1     & 2      & 1     & 3      & ] \\
\end{array}
$
\end{tabular}
\caption{Le format de stockage COO}
\label{fig:coo}
\end{center}
\end{figure}

Dans le format \textit{COO}, les éléments non nuls d'une matrice creuse peuvent être stockés ligne par ligne (voir
l'exemple de la figure~\ref{fig:coo}), colonne par colonne ou dans un ordre quelconque. Cependant, ce format
n'est efficace que pour les matrices creuses non structurées, autrement, nous pouvons remarquer que dans le
cas de stockage ligne par ligne (ou colonne par colonne), le tableau $Lig$ (respectivement, le tableau $Col$)
peut contenir des informations redondantes successives qui pointent le même indice de ligne (respectivement, le même indice de
colonne). Cela impliquerait un espace mémoire non négligeable pour stocker une grande matrice creuse structurée.

\begin{algorithm}
  \SetLine
  \linesnumbered
  \Entree{$Val$, $Lig$ et $Col$ (matrice), $nnz$ (nombre d'éléments non nuls), $n$ (taille de la matrice), $x$ (vecteur)}
  \Sortie{$y$ (vecteur)}
  \BlankLine
  variables $i$, $j$, $id$\; 
  \Pour{$i=0$ \KwA $n-1$}{
    $y[i] \leftarrow 0$\;
  }
  \BlankLine
  \Pour{$id=0$ \KwA $nnz-1$}{
    $i \leftarrow Lig[id]$\;
    $j \leftarrow Col[id]$\;
    $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
  }
\caption{Multiplication matrice-vecteur avec le format COO}
\label{alg:coo}
\end{algorithm}

\subsection{CSR vs. CSC}
Les formats \textit{CSR} (Compressed Sparse Row) et \textit{CSC} (Compressed Sparse Column) sont, probablement,
les plus connus et les plus utilisés pour le stockage des matrices creuses. En effet, ce
sont deux formats qui ne posent aucune condition sur la structure creuse de la matrice et
ne permettent de stocker que les informations nécessaires sur la matrice.
Ils permettent de représenter la matrice creuse sous forme de trois tableaux de données. Comme
le format \textit{COO}, le format \textit{CSR} (respectivement, \textit{CSC}) stocke les $nnz$ éléments non nuls de la matrice
et leurs indices de colonnes (respectivement, leurs indices de lignes) ligne par ligne (respectivement, colonne
par colonne) dans deux tableaux distincts $Val$ et $Col$ (respectivement, $Val$ et $Lig$). Un troisième
tableau d'une taille $n+1$, $Ptr$, est utilisé pour stocker la position du premier élément
non nul de chaque ligne de la matrice (respectivement, de chaque colonne de la matrice). La figure~\ref{fig:csr}
donne les formats de stockage \textit{CSR} et \textit{CSC} de la matrice creuse $A$ donnée comme exemple.
Les algorithmes~\ref{alg:csr} et~\ref{alg:csc} présentent une multiplication
matrice-vecteur avec le format \textit{CSR} et celle avec le format \textit{CSC}, respectivement.

\begin{figure}[!h]
\begin{center}
\begin{tabular}{c}
$A=\left[
\begin{matrix}
a_{00} & 0      & 0     & a_{03} \\
0     & 0      & a_{12} & 0 \\
a_{20} & a_{21} & a_{22} & 0 \\
0     & a_{31} & 0      & a_{33} \\
\end{matrix}
\right]
$
\\
\\
$
CSR:\left\{
\begin{array}{ccccccccccccc}
Val & = & [ & a_{00} & a_{03} & a_{12} & a_{20} & a_{21} & a_{22} & a_{31} & a_{33} &   & ] \\
Col & = & [ & 0     & 3      & 2     & 0      & 1     & 2      & 1     & 3      &   & ] \\
Ptr & = & [ & 0     &        & 2     & 3      &       &        & 6     &        & 8 & ] \\ 
\end{array}
\right.
$
\\
\\
$
CSC:\left\{
\begin{array}{ccccccccccccc}
Val & = & [ & a_{00} & a_{20} & a_{21} & a_{31} & a_{12} & a_{22} & a_{03} & a_{33} &   & ] \\
Lig & = & [ & 0     & 2      & 2     & 3      & 1     & 2      & 0     & 3      &   & ] \\
Ptr & = & [ & 0     &        & 2     &        & 4     &        & 6     &        & 8 & ] \\ 
\end{array}
\right.
$
\end{tabular}
\caption{Le format de stockage CSR et CSC}
\label{fig:csr}
\end{center}
\end{figure}

Les formats \textit{CSR} et \textit{CSC} sont sûrement les plus économes en espace mémoire pour
le stockage des matrices creuses mais, ils ne sont pas forcément les plus efficaces.
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
multiplication matrice-vecteur \textit{CSR} nécessite le chargement de l'indice de colonne
et un accès indirect à la mémoire pour lire la valeur du vecteur $x$ correspondante.
De même, pour le format \textit{CSC} (algorithme~\ref{alg:csc}), chaque opération atomique
de la multiplication requiert un chargement de l'indice de ligne et puis un accès
indirect à la mémoire pour écrire la valeur du vecteur $y$ à l'adresse
correspondante. Donc, la multiplication sur des matrices ayant des lignes ou des
colonnes très creuses engendrera des accès mémoires très irréguliers, ce qui dégaradera
certainement sa performance de calcul.

\begin{algorithm}
  \SetLine
  \linesnumbered
  \Entree{$Val$, $Col$ et $Ptr$ (matrice), $n$ (taille de la matrice), $x$ (vecteur)}
  \Sortie{$y$ (vecteur)}
  \BlankLine
   variables $i$, $j$, $id$\; 
  \Pour{$i=0$ \KwA $n-1$}{
    $y[i] \leftarrow 0$\;
  }
  \BlankLine
  \Pour{$i=0$ \KwA $n-1$}{
    \Pour{$id=Ptr[i]$ \KwA $Ptr[i+1]-1$}{
      $j \leftarrow Col[id]$\;
      $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
    }
  }
\caption{Multiplication matrice-vecteur avec le format CSR}
\label{alg:csr}
\end{algorithm}

\begin{algorithm}
  \SetLine
  \linesnumbered
  \Entree{$Val$, $Lig$ et $Ptr$ (matrice), $n$ (taille de la matrice), $x$ (vecteur)}
  \Sortie{$y$ (vecteur)}
  \BlankLine
  variables $i$, $j$, $id$\; 
  \Pour{$i=0$ \KwA $n-1$}{
    $y[i] \leftarrow 0$\;
  }
  \BlankLine
  \Pour{$j=0$ \KwA $n-1$}{
    \Pour{$id=Ptr[j]$ \KwA $Ptr[j+1]-1$}{
      $i \leftarrow Lig[id]$\;
      $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
    }
  }
\caption{Multiplication matrice-vecteur avec le format CSC}
\label{alg:csc}
\end{algorithm}

\subsection{ELLPACK/ITPACK}
\textit{ELLPACK/ITPACK}, ou plus couramment \textit{ELL}, est parmi les structures de données
les plus génériques pour le stockage des matrices creuses dans la mémoire. De plus, elle est plus
particulièrement adaptée à la représentation des matrices, plus ou moins, structurées et aux
calculs sur des architectures vectorielles. Le stockage en format \textit{ELL} d'une matrice creuse,
$A$, de taille $(n\times n)$, ayant au maximum $m$ éléments non nuls par ligne, utilise deux matrices
de données de taille $(n\times m)$: une matrice de réels, $Val$, et une matrice d'entiers, $Col$.
Une ligne d'indice $i$ de chacune des matrices $Val$ et $Col$ correspond aux éléments non nuls et
leurs indices de colonne, respectivement, de la ligne d'indice $i$ de la matrice creuse $A$.
Les lignes qui ont moins de $m$ éléments non nuls dans la matrice creuse $A$ sont complétées par
des zéros dans les matrices $Val$ et $Col$. La figure~\ref{fig:ell} montre un stockage \textit{ELL}
d'une matrice creuse de taille $(4\times 4)$. Dans cet exemple, les matrices $Val$ et $Col$ sont
complétées par des étoiles ($*$), qui correspondent aux zéros dans la mémoire. L'algorithme~\ref{alg:ell}
montre la mutiplication matrice-vecteur
avec le format \textit{ELL}, telles que les deux matrices $Val$ et $Col$ soient stockées ligne par ligne
dans la mémoire. 
\begin{figure}[!h]
\begin{center}
\begin{tabular}{c}
\begin{tabular}{c}
$A=\left[
\begin{matrix}
a_{00} & 0      & 0     & a_{03} \\
0     & 0      & a_{12} & 0 \\
a_{20} & a_{21} & a_{22} & 0 \\
0     & a_{31} & 0      & a_{33} \\
\end{matrix}
\right]
$
\end{tabular}
\\
\\
\begin{tabular}{cc}
$Val=\left[
\begin{matrix}
a_{00} & a_{03} & * \\
a_{12} & *      & * \\
a_{20} & a_{21} & a_{22} \\
a_{31} & a_{33} & * \\
\end{matrix}
\right],
$
&
$Col=\left[
\begin{matrix}
0 & 3 & * \\
2 & * & * \\
0 & 1 & 2 \\
1 & 3 & * \\
\end{matrix}
\right]
$
\end{tabular}
\end{tabular}
\caption{Le format de stockage ELL}
\label{fig:ell}
\end{center}
\end{figure}

\begin{algorithm}[H]
  \SetLine
  \linesnumbered
  \Entree{$Val$ et $Col$ (matrice), $n$ (taille de la matrice), $m$ (nombre maximum d'éléments non nuls par ligne), $x$ (vecteur)}
  \Sortie{$y$ (vecteur)}
  \BlankLine
   variables $i$, $j$, $l$, $id$\; 
  \Pour{$i=0$ \KwA $n-1$}{
    $y[i] \leftarrow 0$\;
  }
  \BlankLine
  \Pour{$i=0$ \KwA $n-1$}{
    \Pour{$l=0$ \KwA $m-1$}{
      $id \leftarrow n \times i + l$\;
      $j \leftarrow Col[id]$\;
      $y[i] \leftarrow y[i] + Val[id] \times x[j]$\;
    }
  }
\caption{Multiplication matrice-vecteur avec le format ELL}
\label{alg:ell}
\end{algorithm}

\subsection{HYB}
Le format hybride (\textit{HYB}) est la combinaison de deux structures de données, qui sont le format 
\textit{ELL} et le format \textit{COO}. Bien que le format de stockage \textit{ELL} soint bien adapté aux
multiplications matrice-vecteur calculées sur des architectures vectorielles, son efficacité se dégrade
rapidement lorsque le nombre d'éléments non nuls varie substantiellemnt d'une ligne de matrice à une
autre, ce qui est souvent le cas avec les matrices creuses non structurées. En revanche, la performance
des multiplications matrice-vecteur basées sur le format de stockage \textit{COO} est insensible aux
variations du nombre d'éléments non nuls par ligne de la matrice creuse.

Le principe de la structure de données \textit{HYB} consiste à stocker un nombre typique $m$ d'éléments
non nuls par ligne en format \textit{ELL} et le reste des éléments non nuls des lignes exeptionnelles
en format \textit{COO}. Le nombre typique $m$ est souvent déterminé a priori de telle sorte que la
majorité de la matrice creuse ait une structure régulière et ainsi, soit facile à représenter en format
\textit{ELL}. De cette manière, le format \textit{HYB} hérite l'efficacité du format \textit{ELL} due
à ses accès réguliers à la mémoire et la flexibilité du format \textit{COO} qui est insensible à la
structure creuse de la matrice. La figure~\ref{fig:hyb} illustre la structure de données \textit{HYB}
d'une matrice creuse $A$ d'ordre $(4\times 4)$, utilisant un format \textit{ELL} avec 2 éléments par
ligne ($m=2$). L'algorithme~\ref{alg:hyb} présente le code d'une
multiplication matrice-vecteur avec le format de stockage \textit{HYB}.

\begin{figure}[!h]
\begin{center}
\begin{tabular}{c}

\begin{tabular}{c}
$A=\left[
\begin{matrix}
a_{00} & 0      & 0     & a_{03} \\
0     & 0      & a_{12} & 0 \\
a_{20} & a_{21} & a_{22} & 0 \\
0     & a_{31} & 0      & a_{33} \\
\end{matrix}
\right]
$
\end{tabular}
\\
\\
\begin{tabular}{cc}
$
ELL:
$
&
\begin{tabular}{ccc}
$
Val_{ell}=\left[
\begin{matrix}
a_{00} & a_{03}\\
a_{12} & *    \\
a_{20} & a_{21} \\
a_{31} & a_{33} \\
\end{matrix}
\right],
$
&
$
Col_{ell}=\left[
\begin{matrix}
0 & 3\\
2 & * \\
0 & 1 \\
1 & 3 \\
\end{matrix}
\right]
$
\end{tabular}
\\
\\
$
COO:
$
&
\begin{tabular}{ccc}
$Val_{coo}=[a_{22}],$ & $Lig_{coo}=[2],$ & $Col_{coo} = [2]$\\
\end{tabular}
\\
\end{tabular}
\end{tabular}
\caption{Le format de stockage HYB}
\label{fig:hyb}
\end{center}
\end{figure}

\begin{algorithm}[H]
  \SetLine
  \linesnumbered
  \Entree{$Val_ {ell}$ et $Col_{ell}$ (matrice ELL), $Val_ {coo}$, $Lig_{coo}$ et $Col_{coo}$ (matrice COO), 
$n$ (taille de la matrice ELL), $m$ (nombre maximum d'éléments non nuls par ligne dans la matrice ELL), 
$nnz$ (nombre d'éléments non nuls dans la matrice COO), $x$ (vecteur)}
  \Sortie{$y$ (vecteur)}
  \BlankLine
  variables $i$, $j$, $l$, $id$\; 
  \Pour{$i=0$ \KwA $n-1$}{
    $y[i] \leftarrow 0$\;
  }
  \BlankLine
  \Pour{$i=0$ \KwA $n-1$}{
    \Pour{$l=0$ \KwA $m-1$}{
      $id \leftarrow n \times i + l$\;
      $j \leftarrow Col_{ell}[id]$\;
      $y[i] \leftarrow y[i] + Val_{ell}[id] \times x[j]$\;
    }
  }
  \BlankLine
  \Pour{$id=0$ \KwA $nnz-1$}{
    $i \leftarrow Lig_{coo}[id]$\;
    $j \leftarrow Col_{coo}[id]$\;
    $y[i] \leftarrow y[i] + Val_{coo}[id] \times x[j]$\;
  }
\caption{Multiplication matrice-vecteur avec le format HYB}
\label{alg:hyb}
\end{algorithm}
%%-------------------------------------------------------------------------------------------------------%%

%%-------------------------------------------------------------------------------------------------------%%
\section{Parallélisation des méthodes itératives}
\label{sec:parallelisation}
L'adaptation des solveurs itératifs aux architectures des calculateurs parallèles est devenue une
nécessité impérieuse pour la résolution des systèmes linéaires creux de taille toujours croissante.
En effet, les calculateurs parallèles fournissent les ressources et la puissance de calcul nécessaires
pour une résolution haute performance de grands systèmes linéaires creux. 

La parallélisation d'un solveur itératif sur une plateforme parallèle, comportant $p$ processeurs, requiert
tout d'abord un partitionnement des données du système linéaire à résoudre sur l'ensemble des $p$ processeurs.
Le partitionnement de données consiste à attribuer aux différents processeurs des portions, plus ou moins égales,
des vecteurs et des matrices impliqués dans le solveur itératif. Par la suite, tous les processeurs procèdent,
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
itérative mais sur des données différentes. Toutefois, les calculs locaux des différents processeurs sont liés
par des dépendances de données, qui permettent la résolution du système global. A cet effet, des points
de synchronisation doivent être mis en place dans un solveur itératif parallèle, lors desquels des données partagées
sont échangées entre des processeurs voisins.

Dans ce document, nous nous intéressons à la parallélisation des méthodes itératives sur des plateformes de calcul
parallèle à mémoire distribuée. Donc, lors des points de synchronisation, les dépendances de
données sont gérées par des communications à passage de messages entre les processeurs. Les algorithmes
parallèles des méthodes itératives peuvent être classifiés selon la nature \textit{synchrone} ou \textit{asynchrone}
de leurs itérations ainsi que celle de leurs communications~\cite{ref18}.

\subsection{Méthodes itératives SISC}
Les solveurs itératifs parallèles les plus classiques sont ceux dits à Itérations Synchrones
et Communications Synchrones (\textit{SISC}). Le synchronisme des itérations décline du fait que chaque
processeur ne peut commencer le calcul de sa nouvelle itération que lorsqu'il reçoit, de la
part de tous ses voisins, les données partagées calculées à l'itération précédente. Ainsi, tous
les processeurs commencent le calcul de la même itération au même temps et effectuent les échanges
de données partagées à la fin de chaque itération par le biais des communications globales synchrones.

\begin{figure}[!h]
  \centering
  \includegraphics[width=170mm,keepaspectratio]{Figures/SISC}
  \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle SISC avec deux processeurs}
  \label{fig:SISC}
\end{figure}

La figure~\ref{fig:SISC} montre un schéma d'exécution d'un algorithme parallèle \textit{SISC}.
Le synchronisme des itérations des solveurs itératifs parallèles \textit{SISC} impliquent,
exactement, le même nombre d'itérations que leurs homologues en séquentiel. Par conséquent,
les conditions de leur convergence sont faciles à déterminer, vu qu'elles sont identiques
à celles des solveurs itératifs séquentiels. Cependant, l'utilisation des communications
synchrones peut pénaliser les performances de ce type d'algorithmes. En fait, les
communications synchrones engendrent, souvent, des temps d'inactivité des processeurs
(les espaces blancs entre deux itérations calculées par un même processeur dans la
figure~\ref{fig:SISC}), dues soit à l'attente pour un processeur (ou des processeurs)
qu'il soit prêt pour communiquer ou à la vitesse lente des communications synchrones
elles-mêmes dans la plateforme de calcul parallèle.

\subsection{Méthodes itératives SIAC}
Un autre type de solveurs itératifs parallèles a été conçu pour améliorer les performances
de résolution sur les plateformes à réseau d'interconnection lent et/ou hétérogène. Ce sont
les solveurs à Itérations Synchrones et Communications Asynchrones
(\textit{SIAC}). Il permettent de réduire les temps d'inactivité des processeurs entre deux itérations
successives, en utilisant des communications asynchrones et cela, tout en maintenant le
principe des itérations synchrones. En effet, comme les algorithmes \textit{SISC}, un processeur
attend toujours la réception de toutes les données partagées, calculées par ses voisins à
l'itération précédente, avant de commencer les calculs de la nouvelle itération. Toutefois,
les communications synchrones globales utilisées, pour les échanges de données partagées, dans
les algorithmes \textit{SISC} sont remplacées par des envois asynchrones et des réceptions bloquantes.

\begin{figure}[!h]
  \centering
  \includegraphics[width=170mm,keepaspectratio]{Figures/SIAC}
  \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle SIAC avec deux processeurs}
  \label{fig:SIAC}
\end{figure}

La figure~\ref{fig:SIAC} montre un exemple de schéma d'exécution d'un solveur itératif parallèle
\textit{SIAC}. Les conditions de convergence des solveurs itératifs \textit{SIAC} sont identiques à ceux des solveurs
itératifs \textit{SISC} et ainsi, à ceux des solveurs itératifs séquentiels. Cependant, les communications
asynchrones permettent d'effectuer des chauvechaments entre les calculs et les échanges de données.
En fait, un processeur peut envoyer les données partagées à son voisin dès qu'elles soient 
prêtes à être utilisées dans les calculs de la nouvelle itération, ce qui permet de réduire les temps
d'attente de réception de données entre deux itérations successives.

\subsection{Méthodes itératives AIAC}
Les méthodes itératives parallèles synchrones de type \textit{SISC} ou \textit{SIAC} ont
été, largement, utilisées pour la résolution de grands systèmes linéaires creux. Elles sont, relativement,
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
de leurs itérations et/ou de leurs communications pénalise sévèrement les performances de
résolution sur des grilles de calcul à clusters géographiquement distants. En effet, ce
type d'architectures parallèles est souvent caractérisé par la latence forte de ses liens de
communications (plus précisément les liens inter-clusters) et par l'hétérogénéité de ses
ressources. 

Les solveurs itératifs parallèles conçus pour améliorer les performances de résolution
sur les plateformes à ressources hétérogènes et géographiquement distantes sont ceux
dits à Itérations Asynchrones et Communications Asynchrones (\textit{AIAC}). Dans ce
type de solveurs parallèles, chaque processeur exécute ses propres itérations sans
prendre en compte la progression de l'exécution de celles des autres processeurs. En
effet, un processeur passe d'une itération à l'autre sans attendre l'arrivée des nouvelles
données partagées mises à jour par ses voisins. Il utilise les données locales et les
dernières versions des données partagées disponibles au début du calcul de chaque itération.
De plus, pour assurer l'asynchronisme des itérations, les solveurs itératifs \textit{AIAC}
utilisent des communications asynchrones de type envois et réceptions non bloquants.

\begin{figure}[!h]
  \centering
  \includegraphics[width=170mm,keepaspectratio]{Figures/AIAC}
  \caption{Un exemple de schéma d'exécution d'un solveur itératif parallèle AIAC avec deux processeurs}
  \label{fig:AIAC}
\end{figure}

La figure~\ref{fig:AIAC} présente un schéma d'exécution d'un solveur itératif
parallèle de type \textit{AIAC}. Nous pouvons remarquer que l'asynchronisme des
itérations et des communications permet aux processeurs d'éviter les temps
d'inactivité et d'exécuter des itérations différentes à un moment donné de la
résolution. De ce fait, certains processeurs peuvent être plus rapides dans leurs
calculs et effectuer plus d'itérations que les autres. Cependant, les solveurs
itératifs asynchrones \textit{AIAC} ont des conditions de convergence plus strictes
et effectuent plus d'itérations que les solveurs synchrones \textit{SISC} et \textit{SIAC}.
Par conséquent, ils nécessitent une analyse plus élaborée pour déterminer les
bons indicateurs de convergence.
%%-------------------------------------------------------------------------------------------------------%%

%%-------------------------------------------------------------------------------------------------------%%
\section{Conclusion}
Dans ce chapitre, nous avons commencé par présenter les méthodes de résolution des systèmes linéaires creux.
Premièrement, nous avons abordé l'analyse des méthodes directes qui permettent de trouver la solution exacte en un
nombre d'opérations élémentaires fini. Cependant, pour la résolution des systèmes linéaires creux, nous avons mis
en évidence l'existence d'une opération de remplissage à l'étape de factorisation des matrices creuses. L'opération de
remplissage est très coûteuse en termes de temps d'exécution et d'espace mémoire et, plus précisément, pour
les matrices creuses de grande taille. Ensuite, nous avons décrit les méthodes itératives qui permettent de 
calculer une solution approximative en exécutant des itérations successives d'un même bloc d'opérations
élémentaires. Nous avons cité trois grandes familles de méthodes itératives à savoir: les méthodes 
stationnaires, les méthodes non stationnaires et les méthodes multigrilles. Celles-ci sont plus adaptées à la
résolution des systèmes linéaires de grande taille et plus faciles à paralléliser. 

Puis, nous avons étudié différentes structures de données des matrices creuses en mémoire. Elles permettent
un stockage optimisé des éléments non nuls des matrices creuses, de façon à ce que les accès en lecture/écriture
à la mémoire soient faciles à effectuer. Par conséquence, elles permettent d'améliorer les performances de
calcul des méthodes de résolution. Enfin, nous avons donné les principaux points clés de la parallélisation
des méthodes itératives sur les plateformes de calcul parallèle à mémoire distribuée, afin de résoudre des
systèmes linéaires creux à très grande échelle. Nous avons présenté trois différentes méthodes itératives
parallèles, classifiées selon la nature synchrone ou asynchrone de leurs itérations et/ou de leurs communications.
Les méthodes synchrones sont plus efficaces sur des petits clusters homogènes ayant des liens de communications
rapides, alors que les méthodes asynchrones offrent de meilleures performances sur les plateformes à clusters
hétérogènes et géographiquement distants.

Dans la suite de ce document, nous nous intéresserons à la parallélisation des méthodes itératives sur des
clusters équipés de cartes graphiques GPUs. Nous décrirons les principaux points clés des algorithmes parallèles
de quelques méthodes itératives, les plus utilisées, pour la résolution des systèmes linéaires ou non linéaires
creux. Ensuite, nous effectuerons une comparaison de performance entre les solveurs itératifs parallèles
mis en \oe uvre sur les clusters de GPUs et les mêmes solveurs mis en \oe uvre sur les clusters CPUs.      
%%-------------------------------------------------------------------------------------------------------%%

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "these"
%%% End: