• Accès rapide
  • Page d'accueil
  • International (en)
  • Blog de François
  • Blog de Frédéric
  • Maths
  • Informatique
  • Programmation
  • Nos jeux
  • Jeux vidéo
• Pages spéciales
  • Connexion
  • Forum
  • Auteurs
  • Plan du site
• Répertoires parents
  • Informatique
  • Conception et développement d'une interface C/C++ de solveurs numériques en Calcul Parallèle
• Répertoire courant
  • Nos smileys
  • Conception et développement d'une interface C/C++ de solveurs numériques en Calcul Parallèle
  • Développement d'un solveur de programmes linéaires
  • Calculatrice Scientifique
  • Éditeur de schémas et graphiques en SVG
  • TIPE 2006-2007
  • Nos Programmes
  • Programmation

Conception et développement d'une interface
C/C++ de solveurs numériques
en Calcul Parallèle

Table des matières

Résumé du stage

Au cours de ce stage, je devais concevoir et développer une interface C/C++ de solveurs numériques destinés à la résolution de grands systèmes linéaires issus du domaine de l'aérodynamique. J'ai conçu une hiérarchie de classes permettant la manipulation de matrices creuses ainsi qu'une interface commune pour différents solveurs utilisés par l'équipe SAGE. J'ai enfin réalisé une implémentation que j'ai testée sur les grilles de calculs GRID'5000 [8] de l'IRISA.

Accueil et encadrement

J'ai effectué mon stage au sein de l'Institut de Recherche en Informatique et Systèmes Aléatoires, un centre de recherche INRIA situé à Rennes. J'ai été accueilli par l'équipe de recherche SAGE, dirigée par Jocelyne Erhel. SAGE signifie "Simulations et algorithmes sur des grilles de calcul appliqués à l'environnement", les principaux thèmes abordés étant :

• algorithmes numériques, pour la plupart, algèbre linéaire des grands systèmes creux
• calcul parallèle et calcul sur grille
• applications liées à l'environnement, principalement en hydrogéologie et en géophysique.

J'avais pour maître de stage l'ingénieur de recherche Nadir Soualem et pour co-encadrant le doctorant Désiré Nuentsa.

Remerciements

Je tiens tout d'abord à remercier Jocelyne Erhel qui m'a accepté comme stagiaire dans son équipe. Je suis reconnaissant envers tous les membres de l'équipe SAGE pour leur accueil durant ces trois mois. Mes remerciements vont ensuite à Désiré Nuentsa qui m'a fourni son interface C et apporté son aide lorsque j'en avais besoin. Enfin, je souhaite exprimer ma gratitude envers Nadir Soualem, qui m'a accompagné tout le long de ce stage et le remercie pour la relecture de mon rapport.

1. Introduction

Les grands systèmes linéaires creux interviennent dans de nombreux domaines. L'usage de l'informatique et du calcul parallèle est fondamental pour les résoudre efficacement.

Les matrices utilisées dans ces systèmes sont représentées par des formats spécifiques qui tirent profit de la structure creuse. Dans un premier temps, je présenterai mon travail pour concevoir une hiérarchie de classe sparse_interface<T>. Celle-ci permet la manipulation des différents formats de matrices creuses. Je décrirai divers algorithmes pour les opérations sur ces matrices.

Dans une seconde partie, je m'intéresserai à une hiérarchie de classes permettant de disposer d'une interface commune aux solveurs. Ces classes prendront en argument les matrices sparse_interface<T> décrivant le système. Les points communs à chaque solveur seront regroupés autant que possible. Une attention particulière sera donnée aux solveurs parallèles pour lesquels une utilisation du langage MPI [12] et une expérimentation sur grille de calcul est indispensable.

2. Contexte

De nombreux phénomènes physiques sont régis par ce que l'on appelle des équations aux dérivées partielles (EDP). Trouver une expression analytique de leurs solutions est généralement non trivial et elles sont donc résolues numériquement. Pour cela, on réalise une discrétisation, qui va amener à résoudre des systèmes linéaires d'autant plus grands que la précision demandée est fine. Les matrices décrivant ces systèmes sont généralement creuses, c'est-à-dire qu'elles n'ont qu'un petit pourcentage d'éléments non nuls.

Plusieurs solveurs existent et sont utilisés par l'équipe SAGE. Un travail pour trouver une interface commune avait déjà été commencé :

• Dans la plate-forme H2OLAB [9] développé par l'équipe SAGE.
• Par Désiré Nuentsa, dans son travail pour comparer les performances de divers solveurs.

Il m'a alors été demandé de concevoir une hiérachie de classes étendant le travail déjà effectué, pour obtenir une interface uniforme.

Je devais travailler dans un environnement Linux et développer en C++. Il m'a été demandé de respecter quelques conventions dans la syntaxe du code. Je devais par exemple le documenter à l'aide de Doxygen [4]. Pour les outils de développement, j'ai choisi de coder avec emacs [7] et ai manipulé des outils comme gcc [5], gdb [6] ou valgrind [22]. Il m'a aussi été suggéré d'utiliser CMake [3] pour la gestion de la compilation. Pour partager mon travail, je me suis servi de la forge collaborative de l'INRIA et notamment du gestionnaire de version svn. Le travail de conception a été réalisé à l'aide du logiciel Umbrello [21]. Enfin, MPICH2 [13] a servi d'implémentation pour les calculs en parallèle.

Pour valider mon travail, j'ai réalisé divers programmes de tests. J'ai aussi mis au point des tests de non régression. Enfin, j'ai expérimenté mon programme sur la grille GRID' 5000 [8] de l'IRISA et comparé l'adéquation des résultats avec les implémentations déjà utilisées par l'équipe SAGE.

3. Hiérarchie de classes pour les matrices creuses

3.1. Implémentation de base

3.1.1. Description du travail demandé

Dans un premier temps, je devais concevoir et implémenter une hiérarchie de classes pour la manipulation des matrices creuses. Ces matrices sont carrées de taille size et avec nnz éléments non nuls. Trois formats différents m'ont été présentés :

• Le format COO. Il est formé de trois tableaux de taille nnz : rows, columns et values. Tout k $\left(0\le k<\mathrm{size}\right)$ fournit un élement de coordonnées (rows[k], columns[k]) de valeur values[k].
• Le format CSR (Compressed Sparse Row). Il est formé d'un tableau pointer de taille size+1 et de deux tableaux columns et values de taille nnz. columns et values sont identiques à COO. Le numéro de ligne de l'élément est déterminé grâce au tableau pointer, qui donne la position k du premier élément non nul de la ligne i.
• Le format CSC (Compressed Sparse Column). Idem que pour CSR, mais la matrice se lit en colonne d'abord.

Considérons à titre d'exemple la matrice suivante :

$A=\left(\begin{array}{ccccc}1& 1& -3& 0& 0\\ -2& 5& 0& 0& 0\\ 0& 0& 4& 6& 4\\ -4& 0& 2& 7& 0\\ 0& 8& 0& 0& -5\end{array}\right)$

Les tableaux décrivants la matrice dans ces différents formats sont

• COO :

  rows    = [ 0  0  0  1  1  2  2  2  3  3  3  4  4]
  columns = [ 0  1  2  0  1  2  3  4  0  2  3  1  4]
  values  = [ 1  1 -3 -2  5  4  6  4 -4  2  7  8 -5]    

• CSR :

  pointer = [ 0        3     5        8       11    13]
  columns = [ 0  1  2  0  1  2  3  4  0  2  3  1  4]
  values  = [ 1  1 -3 -2  5  4  6  4 -4  2  7  8 -5]    

• CSC :

  pointer = [ 0        3     5        8       11    13]
  columns = [ 0  1  3  0  1  4  0  2  3  2  3  2  4]
  values  = [ 1 -2 -4  1  5  8 -3  4  2  6  7  4 -5]

Il m'a aussi été demandé de prendre en compte le cas des matrices symétriques. Les structures précédentes peuvent être reprises, mais seule la partie supérieure doit être stockée, ce qui divise presque par deux le nombre d'éléments à mémoriser.

Ces formats devaient être implémentés dans une interface commune. Les matrices devaient notamment posséder des constructeurs permettant de passer d'un format à un autre. En outre il devait être possible d'effectuer différentes opérations : écrire/lire un élément, addition/soustraction de matrices, multiplication matrice-vecteur, multiplication/division d'une matrice par un scalaire, transposée. Une surcharge d'opérateur () pour utiliser des expressions A(i,j) a été suggérée pour plus tard.

Enfin, il a été envisagé une classe permettant de manipuler de manière générique chacun de ces types de matrices creuses.

3.1.2. Première proposition de hiérarchie de classes

Le diagramme de la figure 1 est une première proposition pour la hiérarchie de classes. Elle sera révisée plus loin pour prendre en compte les résultats de l'implémentation. Toutes les classes possèdent un paramètre T qui est le type des coefficients des matrices.

Étant donné que les trois formats de matrices possèdent les mêmes opérations, nous regroupons ces dernières dans une interface UML sparse_interface<T>. Comme une telle notion n'existe pas en C++, celle-ci sera représentée lors de l'implémentation par une classe parente. Remarquons au passage que dans le diagramme, certaines fonctions membres prennent une matrice A en paramètre : ce sont celles qui seront implémentées en friend.

Nous introduisons ensuite trois classes CSC, COO, CSR qui implémentent l'interface sparse_interface<T> et représentent les structures de matrices creuses discutées plus haut. Ces classes possèdent en outre des constructeurs à partir des autres types de façon à permettre les conversions.

Les matrices symétriques sont représentées par des classes filles supplémentaires. Comme les données stockées en CSC et CSR sont exactement les mêmes dans le cas des matrices symétriques, nous nous contentons d'ajouter deux classes COOsym et CSRsym héritant respectivement des classes COO et CSR.

Enfin, nous ajoutons une classe sparse<T> qui est une composition des trois structures précédentes et implémente elle-même l'interface sparse_interface<T>. Une seule de ces structures est activée à un moment donné. Les opérations effectuées sur la classe sparse<T> sont alors répercutées sur la structure active.

3.1.3. Première implémentation

Comme demandé l'implémentation est réalisée en C++ avec utilisation de patrons de classe pour avoir un type de coefficients générique dans les matrices. Néanmoins dans un premier temps j'ai codé en utilisant pour paramètre T = double de façon à éviter les difficultés liées à la syntaxe des patrons. Cette restriction initiale ne changeant pas grand chose puisque la réflexion algorithmique porte principalement sur la structure des matrices. Cette première implémentation ne prend pas encore en compte les classes sparse<T> ou symétriques.

Les coordonnées et valeurs décrivant les matrices creuses sont stockées dans une structure de tableau, ce qui permet d'accéder très rapidement à un élément lors des opérations algébriques ainsi que d'occuper le moins de mémoire possible. Par contre, contrairement à une structure de type liste par exemple, l'insertion de nouveaux éléments risque d'être lent comme nous allons le voir par la suite. Pour représenter ces tableaux, j'ai choisi d'utiliser le type vector de la librairie standard C++. Il permet en effet de conserver une disposition simple en mémoire tout en facilitant l'allocation/désallocation. La documentation avertit effectivement de l'inefficacité de l'opération insert().

Les tableaux utiliss pour chaque type de matrice sont ceux donnés dans la description du travail demandé. A partir de maintenant, nous utiliserons donc nnz et size pour désigner le nombre d'éléments non nuls et la taille de la matrice. Nous utiliserons une variable i ou j pour parcourir des tableaux de taille size+1 (les pointer) et une variable k pour les tableaux de taille nnz (columns, rows et values). Le parcours de la matrice se fait alors en O(nnz):

// COO
Pour k de 0 à nnz - 1
 [opérations sur le k-ième élément]
FinPour

// CSR & CSC
Pour i de 0 à size - 1
  Pour k de pointer[i] à pointer[i+1] - 1
    [opérations sur le k-ième élément]
  FinPour
FinPour

Précisons une propriété de l'implémentation : les éléments de la matrice sont triés selon l'ordre lexicographique des coordonnées (i,j). Pour les matrices COO et CSR, on lit "en ligne d'abord" et pour CSR "en colonne d'abord". Cela peut sembler une contrainte supplémentaire pour l'insertion mais une fois encore, cela permet d'accélerer la recherche d'un élément ainsi que le parcours lors des opérations algébriques.

Intéressons-nous maintenant aux opérations sur les matrices. Les algorithmes d'addition/soustraction et de conversion/transposition seront discutés dans des sections ultérieures. Tout d'abord, notons que les algorithmes de CSR/CSC sont symétriques. Les opérations de multiplications/divisions par un scalaire sont très simples. Il suffit de parcourir values et d'appliquer l'opération sur chacun des éléments. L'opération y = Ax est à peine plus difficile. On commence par faire y = 0 puis on parcourt les éléments de la matrice. Pour COO, on le fait selon k en effectuant à chaque fois l'opération y[rows[k]]+=values[k]*x[columns[k]]. Pour CSR, on fait deux boucles imbriquées en (i,k) (voir plus haut) et on effectue y[i]+=values[k]*x[columns[k]]. Pour CSC, on parcourt en (j,k) et on effectue y[rows[k]]+=values[k]*x[j]. Toutes ces opérations se réalisent donc en O(nnz).

Remarquons ensuite que les opérations getValue et setValue sont très similaires : on recherche la position d'un élément (i,j). Si on le trouve, on retourne/modifie sa valeur. Sinon, on s'arrête à la première position selon l'ordre lexicographique. En effet pour getValue il est inutile de chercher plus loin et pour setValue c'est la position où on doit insérer l'élément. Nous introduisons donc une fonction getPosition pour réaliser cet opération. Pour COO, il s'agit d'une boucle pour k allant de 0 à nnz - 1 et pour CSR une boucle allant de pointer[i] à pointer[i+1]. Les tableaux étant triés, nous pouvons utiliser une recherche dichotomique.

La complexité des opérations getValue/setValue varie beaucoup. Dans le meilleur des cas, on trouve directement l'élément et on le modifie. Dans le pire cas de getPosition, on doit parcourir toute la matrice pour COO (toute la ligne pour CSR) donc une complexité linéaire selon le nombre d'éléments parcourus (ou logarithmique si on utilise la dichotomie). Notons tout de suite que pour CSR/CSC il y a moins d'éléments à parcourir. La recherche est donc plus efficace mais l'amélioration par dichotomie est moindre. Remarquons une nouvelle fois l'opération la plus coûteuse ici : l'insertion d'un élément dans le tableau lors d'un setValue, où il faut en effet décaler tous les éléments qui suivent, avec une complexité linéaire. De plus cette opération nécessite potentiellement une réallocation.

Pour vérifier les performances de la première implémentation, j'ai mesuré le temps d'exécution (sauf pour la somme/soustraction matrice-matrice). La matrice utilisée est une matrice bande de taille 10000 de "largeur" 11, c'est-à-dire de l'ordre de 110000 éléments non nuls, et 11 éléments non nuls par ligne/colonne. Toutes les opérations citées sauf le remplissage se réalise en 0s. Pour ce dernier on obtient, on obtient les résultats donnés dans le tableau de la figure 2.

Nous constatons donc que toutes les opérations (produit matrice-vecteur, multiplication par un scalaire...) sont négligeables devant le remplissage et c'est donc ce dernier que nous devons améliorer. Comme annoncé, la dichotomie améliore les performances pour COO mais a un effet négligeable sur CSR/CSC. Une solution consistant à mémoriser les dernières positions à été envisagée. Toutefois cela accélèrerait seulement la recherche et non le décalage des éléments lors de l'insertion. De plus cela ajoute des hypothèses sur la manière dont l'utilisateur remplit la matrice.

La mesure du temps de remplissage pour une taille de 100000 est de l'ordre de 20 minutes ! C'est bien trop, sachant que l'utilisation sur un poste individuel peut atteindre cet ordre de grandeur.

3.1.4. Algorithme d'addition/soustraction de matrices

Prenons le cas de la soustraction de deux matrices COO. Tout d'abord, on ne souhaite pas réutiliser les opération getValue et setValue qui nécessite un parcours pour chaque coefficient. On va plutôt parcourir en parallèle les matrices A et B grâce à deux curseurs kA et kB. La matrice résultante AmoinsB sera parcouru avec un curseur kAmoinsB. On constate tout d'abord que AmoinsB.nnz ≤ A.nnz + B.nnz. Nous commençons donc par agrandir les tableaux de AmoinsB à la taille A.nnz + B.nnz, de façon à disposer d'assez de place. Les curseurs sont ensuite tous initialisés à 0. Puis on compare les coordonnées de l'élément de A en position kA et de celui de B en position kB. Si elles sont égales on soustrait les deux éléments, on stocke le résultat dans AmoinsB à la position kAmoinsB et on incrémente tous les curseurs. Sinon, on ne prend en compte que l'élément avec la plus petite coordonnée et on réalise la soustraction avec zéro comme autre opérande. Seul le curseur de kAmoinsB et de l'élément choisi sont incrémentés. On continue jusqu'à ce que l'on ait parcouru entièrement les matrices A et B. Si on note $E_A$ (respectivement $E_B$ ) l'ensemble des coordonnées (i, j) des éléments non nuls de A (respectivement B). La complexité dans le pire des cas ( $E_A\cap E_B=\varnothing$ ) est O(A.nnz + B.nnz) et dans le meilleur des cas ( $E_A\subseteq E_B$ ou $E_B\subseteq E_A$ ) en O(max(A.nnz, B.nnz)).

L'algorithme se généralise facilement au cas de l'addition. Il peut aussi s'appliquer aux structures CSR/CSC, il suffit simplement de modifier la façon dont les matrices sont parcourues. Pour cela, on commence par faire une boucle avec une variable i parcourant toutes les lignes/colonnes. Puis on répète l'algorithme précédent en faisant cette fois varier kA et kB sur la ligne/colonne i, c'est-à-dire respectivement de A.pointer[i] à A.pointer[i+1] et B.pointer[i] à B.pointer[i+1].

3.1.5. Algorithmes de conversion et transposition

Dans cette section, nous nous penchons sur le cas des conversions entre les différentes structures. En fait, les algorithmes mis en oeuvre pourront aussi être utilisés pour la transposition. Nous allons voir que tous s'effectuent rapidement : en O(nnz) ou O(nnz + size). Commençons par le cas de la conversion CSR ↔ COO. Nous remarquons en effet que les valeurs des tableaux columns et values sont exactement les mêmes.

Les valeurs de rows sont une succession d'indices i rangés dans l'ordre et répété autant de fois qu'il y a d'éléments dans la ligne. Mais ce nombre d'éléments est exactement donné par pointer[i+1]-pointer[i] on a ainsi une conversion directe CSR à COO :

[complexité: 3*nnz = O(nnz)]
Fonction CSR_TO_COO(CSR A, COO B)

  (* Valeurs identiques *)
  B.nnz := A.nnz
  B.size := A.size
    (* !! Recopie en 2*nnz opérations !! *)
  B.columns := A.columns
  B.values := A.values
  
  (* Remplissage du tableau "rows" *)
  Pour i de 0 à A.size - 1
     Pour k de A.pointer[i] à A.pointer[i+1] - 1
        B.rows[k] := i
     FinPour
  FinPour

FinFonction

Pour l'opération inverse, on parcourt les éléments de rows. Dès que l'on trouve un changement de valeur, c'est que l'on est sur une nouvelle ligne. On peut mémoriser l'indice de ces changements pour en déduire le tableau pointer. Il faut cependant faire attention à bien traiter le cas des lignes vides :

[complexité: 3*nnz + size = O(nnz + size)]
Fonction COO_TO_CSR(COO A, CSR B)
  (* Valeurs identiques *)
  B.nnz := A.nnz
  B.size := A.size
    (* !! Recopie en 2*nnz opérations !! *)
  B.columns := A.columns
  B.values := A.values
  
  (* Remplissage du tableau "pointer" *)
  Si A.nnz > 0
    (* Détermination de B.pointer[i] pour les lignes non vides *)
    i := A.rows[0]
    estVide[i] := FAUX
    B.pointer[i] := 0

    Pour k de 0 à A.nnz - 1
      Si A.rows[k] != i
        i := A.rows[k]
        estVide[i] := FAUX
        B.pointer[i] := k
      FinSi
    FinPour

    (* Remplissage de B.pointer[i] pour les lignes vides *)
    B.pointer[A.size] := A.nnz
    Pour i de A.size - 1 à 0
      Si estVide[i]
        B.pointer[i] := B.pointer[i+1]
      FinSi
    FinPour
  FinSi

FinFonction

Étudions maintenant le cas du passage de "en ligne d'abord" à "en colonne d'abord" (ou l'inverse). Dans ce cas, l'ordre de parcours change totalement. Lors de la recopie, on parcourt la matrice source. Pour se repérer dans la matrice destination, on est obligé de calculer un pointer start[j] sur le début de chaque colonne j (que l'on ne possède pas puisqu'au mieux on pointe sur le début des lignes) ainsi que le nombre d'éléments nb[j] déjà rangés dans cette colonne. La nouvelle position de l'élément sera start[j] + nb[j]. Il y a ensuite des modifications pour chaque cas selon les tableaux que l'on cherche. Par exemple si on veut un CSC, son tableau pointer n'est rien d'autre que start. Voici par exemple l'algorithme pour passer de COO à CSC :

[complexité: size + nnz + size + size + size + nnz = 2*nnz + 4*size = O(nnz + size)]
Fonction COO_TO_CSC(COO A, CSC B)
  (* Valeurs identiques *)
  B.size := A.size
  B.nnz := A.nnz

  (* Determination du nombre d'élements par colonne *)
  Pour j de 0 à A.size - 1
    nb[j] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    nb[A.columns[k]] := nb[A.columns[k]] + 1
  FinPour

  (* Remplissage du tableau "pointer" pour les colonnes non vides *)
  k := 0
  Pour j de 0 à A.size - 1
    Si nb[j] > 0
      B.pointer[j] := k
      k := k + nb[j]
    FinPour
  FinPour

  (* Remplissage du tableau "pointer" pour les colonnes vides *)
  B.pointer[A.size] := A.nnz
  Pour j de A.size - 1 à 0
    Si estVide[j]
      B.pointer[j] := B.pointer[j+1]
    FinSi
  FinPour

  (* Remplissage des tableaux "rows" et "values" *)
  (* nb[j] = nombre d'élément déjà stocké dans la colonne j *)
  Pour j de 0 à A.size - 1
    nb[j] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    j := A.columns[k]
    kB := B.pointer[j] + nb[j]
    B.rows[kB] := A.pointer[k]
    nb[j] := nb[j] + 1
  FinPour

FinFonction

L'algorithme précédent s'applique aussi aux conversions restantes et même à l'opération de transposition (cf Annexe A). Dans un premier temps, je n'avais pas implémenté tous ces algorithmes et ai expérimenté deux méthodes alternatives :

• le passage par un objet intermédiaire. Par exemple CSC_TO_CSR = COO_TO_CSR ∘ CSC_TO_COO
• l'utilisation de l'algorithme de remplissage (voir section suivante). On parcourt la matrice source et on construit la matrice destination.

J'ai mesuré le temps d'exécution des opérations pour une matrice de taille 100000 et une largeur de bande 41. Notons qu'une erreur d'une dizaine de millisecondes peut apparaitre pour une même mesure. On constate que l'utilisation du remplissage par setValue coûte très cher et qu'une conversion directe est préférable.

1. Conversion directe COO ↔ CSR et COO ↔ CSC et utilisation de setValue pour les autres.

   COO→CSR	COO → CSC	CSR→COO	CSR → CSC	CSC → COO	CSC → CSR
   60ms	270ms	90ms	8140ms	300ms	8460ms

2. Conversion directe COO ↔ CSR et COO ↔ CSC / passage par COO intermédiaire pour les autres.

   COO → CSR	COO → CSC	CSR → COO	CSR → CSC	CSC → COO	CSC → CSR
   60ms	260ms	100ms	400ms	300ms	400ms

3. Conversion directe.

   COO → CSR	COO → CSC	CSR → COO	CSR → CSC	CSC → COO	CSC → CSR
   80ms	250ms	90ms	260ms	300ms	260ms

Les coût empiriques de ces conversions/transpositions seront analysés plus en détail dans la section 3.2.3.

3.2. Amélioration de l'implémentation

3.2.1. Amélioration du remplissage des matrices

Pour trouver une solution au problème du remplissage, je me suis intéressé à la librairie open source MTL [14] qui implémente les formats CSR/CSC. Dans sa documentation, il est précisé que la modification directe d'un élément n'est pas possible pour des raisons de performance. L'opération de remplissage se décompose grosso modo en trois étapes :

1. Initialisation : on réorganise le tableau pour que chaque ligne soit stockée dans un bloc de taille minimale slot_size.
2. Édition : l'utilisateur apporte des modifications. On écrase directement un élément s'il existe déjà. Lorsqu'il doit être inséré, on effectue l'opération s'il reste de la place dans le bloc de sa ligne. Sinon on mémorise le triple (i, j,valeur).
3. Finalisation : on donne la taille finale pour les blocs de chaque ligne et on effectue les opérations mémorisées.

Par rapport à la méthode naïve, on ne déplace de "gros" blocs qu'à l'étape 1 et 3. De plus ces déplacements peuvent se faire avec un décalage strictement plus grand que 1. Lorsque l'on insère un élément (soit à l'étape 2 s'il y a de la place, soit à l'étape 3 si cet élément a été mémorisé) on ne décale que les éléments à l'intérieur d'un bloc. Ceux-ci correspondant au nombre d'éléments non nuls d'une ligne, on peut espérer qu'ils sont peu nombreux. Dans le pire des cas, ils sont au nombre de size. Le gain est assez impressionnant : pour la matrice vue dans la section 3.1.3 mais avec une taille de 1 million, le remplissage avec la librairie MTL ne met que 12 secondes ! J'ai donc cherché à reprendre cet algorithme et ai obtenu un temps similaire : 10 secondes.

Dans la librairie MTL, l'utilisateur doit manipuler un objet supplémentaire inserter. Les étapes 1) et 3) se font alors respectivement à la création et à la destruction de cet objet. A l'étape 2) l'écriture A.setValue(i, j, value) est remplacée par ins[i][j] << value, où ins est un inserter pour la matrice A. Cela nécessite des surcharges d'opérateurs pour permettre une écriture spéciale et la manipulation d'un objet intermédiaire. Tout ceci ne me semblant pas très intuitif pour l'utilisateur, j'ai opté pour l'ajout de fonctions membres begin() et end() qui annonce le début et la fin de l'édition. Ce genre de technique se retrouve dans d'autres applications informatiques (sémaphores dans les systèmes d'exploitation, annonce de transaction dans les bases de données...) et me semble donc plus naturel. Le remplissage de la matrice A se fera dorénavant de la façon suivante :

// étape 1
A.begin();

// étape 2 = appels à setValue
 ...
A.setValue(i, j, v);
 ...

// étape 3
A.end();

Cette modification pose la question de savoir si on autorise des opérations durant la phase de remplissage i.e. entre les appels à begin() et end(). Bien qu'il serait possible d'ajouter des traitements spécifiques permettant la modification des valeurs mémorisées durant le remplissage, cela alourdirait le code. On pourait penser autoriser des opérations simple comme la multiplication par un scalaire (simple application sur chaque élément) par rapport à la somme (deux matrices à prendre en compte avec des états potentiellement différents). Néanmoins cela risque d'entrainer une confusion pour l'utilisateur : savoir quels objets ont une séquence d'édition ouverte, quelles opérations sont autorisées etc. Par conséquent on préfère n'autoriser que deux opérations entre begin() et end() : getValue et setValue.

Notons aussi que la puissance de l'algorithme réside dans le découpage en blocs de lignes (ou de colonnes). Pour cela, il est nécessaire de disposer du tableau pointer. L'algorithme ne semble donc pas applicable à la structure COO. Néanmoins on peut passer très rapidement de la forme COO à la forme CSR comme nous l'avons vu dans la section précédente. Ainsi effectuer ces conversions aux étapes 1 et 3 apportera un coût supplémentaire négligeable par rapport au reste des opérations.

Passons maintenant à une description plus précise de l'algorithme dans le cas de CSR. Chaque ligne de la matrice se voit allouer un certain espace entre pointer[i] et pointer[i + 1] dans lequel vont être stockés les nouveaux éléments (les valeurs de columns et values). L'étape 1) consiste à agrandir cet espace s'il est inférieur à une valeur slot_size (par défaut à 5). Après l'étape 1, toutes les cases ne sont pas forcément occupés. Ainsi les éléments de la ligne i sont stockés entre pointer[i] et pointerEnd[i] (pointerEnd étant un nouvel objet), les autres étant des cases libres. A l'étape 2) on ajoute des éléments et on risque de dépasser la place réservée. Dans ce cas, on les mémorise dans une structure de la librairie standard, une map toInsert [23]. La structure de la mémoire durant l'édition est représentée sur la figure 3.

Notons que l'implémentation de ce nouvel algorithme ajoute un tableau pointerEnd et une map toInsert. Ces objets sont aloués dynamiquement et libérés à la fin de l'édition. Voici le pseudo-code des fonctions begin, setValue et end. Comme on souhaite simplement donner l'idée, on ne détaille pas toutes les opérations. Les fonctions copy et copy_backward de la librairie standard sont utilisées pour une meilleur performance.

Fonction begin

  (* création des tableux "pointerEnd" et "new_pointer" de taille size + 1 *)
   new_pointer[0] := 0

   (* Détermination des nouvelles tailles/positions des blocs *)
  Pour i de 0 à size - 1
    nombre_elements_dans_ligne_i := pointer[i+1] - pointer[i]
    pointerEnd[i] := new_pointer[i] + nombre_elements_dans_ligne_i
    new_pointer[i+1] := new_pointer[i] + max(nombre_elements_dans_ligne_i, slot_size)
  FinPour

  (* Redimensionnement des tableaux values et columns *)

  (* Déplacement des blocs : comme on ne diminue jamais la taille des blocs, on les
déplace tous vers la fin. On commence donc à déplacer les blocs finaux. *)
  Pour i de size - 1 à 0
     (* Déplacement du bloc i : recopie par la fin, pour éviter les chevauchements *)
  FinPour

   (* Remplacement de pointer *)
   pointer := new_pointer

FinFonction

Fonction setValue(i, j, value)
  
  (* Recherche de la position *)
  (position, trouvé) := getPosition(i, j)

  Si trouvé
     (* Mise à jour de la valeur dans le tableau *)
  Sinon
     Si pointerEnd[i] != pointer[i + 1]
        (* Il reste de la place dans le bloc i *)
        (* Insertion de la valeur en recopiant le reste de la ligne par la fin *)
     Sinon
        (* Stockage dans la map *)

         (position, trouvé) := toInsert.find(i,j)
         Si trouvé
           (* Mise à jour de la valeur dans la map *)
         Sinon
           (* Insertion d'une nouvelle valeur dans la map *)
         FinSi 

     FinSinon
  FinSinon

FinFonction

Fonction end

   (* Création du tableau "new_pointer" de taille size + 1 *)
   new_pointer[0] := 0
   it := Itérateur_sur_le_début(toInsert)

   (* Détermination des nouvelles tailles/positions des blocs *)
   Pour i de 0 à size - 1

     (* Détermination du nombre d'élément dans la ligne i *)
     nombre_elements_dans_ligne_i := pointerEnd[i+1] - pointer[i]

     TantQue "it n'est pas à la fin de la map et que l'élément qu'il pointe est sur la ligne i"
        nombre_elements_dans_ligne_i := nombre_elements_dans_ligne_i + 1
        Avancer(it)
     FinTantQue 

     new_pointer[i+1] = new_pointer[i] + nombre_elements_dans_ligne_i

   FinPour 

   (* Si nnz a augmenté, on redimensionne "columns" et "values" maintenant *)     

   (* Déplacement des blocs *)

   (* Si nnz a diminué, on redimensionne "columns" et "values" maintenant *)     

   (* Détermination du nouveau pointerEnd *)
   Pour i de 0 à size - 1
     pointerEnd[i] = new_pointer[i] + pointerEnd[i] - pointer[i]
   FinPour

   (* Remplacement de pointer *)
   pointer := new_pointer

   (* Insertion des éléments stockées dans toInsert *)
  
FinFonction

Dès la première implémentation, on pouvait constater que les classes COO, CSR et CSC ont beaucoup de code en commun. L'ajout de ce nouvel algorithme de remplissage a accentué encore davantage le phénomène. Ainsi, nous insérons dans notre hiérarchie une classe sparse_common qui implémente sparse_interface<T> et regroupe ce code commun. Les classes COO, CSR et CSC héritent maintenant de sparse_common.

3.2.2. Étude du temps de remplissage en fonction du slot_size

Dans cette section, nous allons observer l'influence du slot_size. Il s'agit de la taille minimale des blocs alloués pour chaque ligne. Il peut être passé en paramètre à la fonction begin() par l'utilisateur. Plaçons nous dans le cas où le nombre d'éléments non nuls par ligne est constant. Si le slot_size est inférieur à cette valeur, on va devoir utiliser la map. Si il est supérieur, on déplace des blocs pour rien et on agrandit trop les tableaux. C'est pourquoi la documentation de MTL propose à l'utilisateur de préciser ce paramètre si possible. Lorsqu'il varie, il propose de fournir "une valeur supérieur à la moyenne du nombre d'élément non nuls par ligne et peut-être inférieur au maximum".

Nous prenons tout d'abord le cas de notre matrice bande de taille 100000 et de largeur de bande 41. Son nombre d'éléments non nuls par ligne est (presque) constant de valeur 41. La courbe de la figure 4 donne le temps de remplissage (en millisecondes) pour un slot_size qui varie. On observe une évolution linéaire avec deux fragments de courbes : une décroissance élevée de 0 à 41 et une croissance faible de 41 à 100. Comme annoncé, le slot_size donnant la performance optimale est celui correspondant au nombre d'éléments non nuls de la matrice. Notons qu'il vaut mieux surestimer le slot_size plutôt que le sous-estimer car la pente est plus faible.

Considérons maintenant une matrice de taille 50000 avec un nombre d'éléments non nuls variant de 50 à 200 (figure 5). La moyenne mesurée du nombre d'éléments non nuls par ligne est 124,975. Nous obtenons cette fois-ci un changement de pente plus doux, entre la moyenne et le maximum. Comme précédemment, la courbe décroit vite avant slot_size = 125 et reste pratiquement minimale au niveau de slot_size = 200 (1280ms). Elle croit ensuite très lentement (1290ms pour slot_size = 210, 1310ms pour slot_size = 225). A nouveau il faut surestimer le slot_size. De plus, il vaut mieux ici être prêt du nombre maximum d'éléments non nuls.

3.2.3. Mesure des performances

Pour mesurer ces performances, nous utilisons une matrice bande de largeur de bande 21. Pour la somme matrice-matrice, nous utilisons pour second opérande une matrice formée d'une bande horizontale de largeur 21. Les valeurs de nnz sont donc environ de 21 * size. Nous étudions le temps en millisecondes mis pour chaque opération en faisant varier la taille de la matrice. Pour une taille donnée, on observe que la mesure peut varier d'une dizaine de millisecondes. Par conséquent chaque mesure a été effectuée en triple pour tenir compte de cette erreur. Les résultats sont donnés en Annexe C (figure 12).

Tout d'abord, nous avons trouvé pour pratiquement tous nos algorithmes une complexité de la forme a*nnz + b*size. Comme ici nnz est pratiquement proportionnel à size, il est normal que l'on observe des droites. Nous remarquons que c'est aussi le cas pour l'algorithme de remplissage (12.4). Ce dernier se distingue bien des autres par son coût : pour la taille maximale, il dure une trentaine de secondes au lieu d'une seconde pour les autres algorithmes. Les trois premiers schémas (12.1, 12.2 et 12.3) comparent ces opérations peu coûteuses pour COO, CSR et CSC. On remarque que l'on a toujours l'ordre suivant pour les coûts : produit matrice-vecteur (nnz opérations), somme matrice-matrice (A.nnz + B.nnz = 2nnz opérations), transposition (2nnz + 4size opérations). Le coût des conversions varie mais il est généralement proche de celui de la transposition ou du produit matrice-vecteur.

Maintenant si on considère une comparaison entre les différents types de structures, nous observons que pour le remplissage (12.4), la transposition (12.7) et le produit matrix-vecteur (12.5) les droites se superposent quasiment. Ainsi pour ces méthodes, nous pouvons indifféremment utiliser COO, CSR, CSC même si la première semble légèrement plus coûteuse. Pour la somme matrice-matrice (12.6), CSR/CSC sont équivalentes mais COO est plus coûteuse. Une explication probable est que, dans l'algorithme de sommation, la comparaison de deux coefficients non nuls des matrices A et B est plus complexe dans ce cas. En effet, en CSR/CSC on est à chaque instant à une ligne/colonne fixée et il suffit de comparer la colonne/ligne. A l'inverse dans COO, on doit à chaque fois réaliser une comparaison lexicographique sur les couples (ligne, colonne).

La principale différence réside donc dans les complexités des conversions (12.8). Remarquons que toutes les conversions autres que COO ↔ CSR ont pratiquement la même pente. En fait, on retrouve la pente de la transposition (12.7) : environ 100ms en 100000 et 1300ms en 1000000. A l'inverse, la pente des conversions COO ↔ CSR est plus faible. Comme nous l'avions vu dans la section 3.1.5, les algorithmes sont plus simples lorsqu'il n'est pas utile de changer l'ordre de lecture de la matrice. Remarquons toutefois un résultat pouvant sembler surprenant. La conversion CSR → COO coûte plus cher que COO → CSR, alors que l'algorithme semble plus simple. En fait, si on observe attentivement, chaque fois que le résultat d'une conversion/transposition est COO, la pente est légèrement plus élevée. Dans les calculs de complexité, nous n'avons pas pris le coût de l'allocation/remplissage des nouveaux tableaux. Or pour COO tous les tableaux sont de taille nnz mais pour CSC/CSR le tableau pointer ne contient que size éléments, ce qui peut expliquer la performance moindre.

3.2.4. Implémentation des matrices symétriques

Nous souhaitons maitenant utiliser des objets particuliers pour les matrices symétriques. Comme annoncé dans la section 3.1.2, il suffit de réaliser deux classes COOsym et CSRsym héritant respectivement de COO et CSR. Les tableaux utilisés sont les mêmes, mais seuls les éléments de la partie supérieure sont stockés.

Il convient donc de vérifier si l'on doit réimplémenter les opérations. Les opérations setValue/getValue sur le coefficients de coordonnées (i, j) fonctionne toujours si $i\le j$ . Dans le cas contraire, il suffit d'appeler l'opération sur le coefficient (j, i). Une manière simple de faire cela est donc de ré-implémenter getValue et setValue en prenant en compte ces deux cas et en appelant les opérations de la classe mère avec (i, j) ou (j, i). Notons au passage que tout se passe comme si le setValue modifiait simultanément deux éléments dans le cas des éléments extra-diagonaux.

Les opérations de multiplication/division par un scalaire ou d'addition/soustraction matrice-matrice n'ont pas besoin d'être ré-implémentées dans les classes filles. On s'assure rapidement que l'ensemble des matrices symétriques de taille n est bien stable pour ces opérations.

Il reste enfin le cas de la multiplication matrice vecteur y = Ax. Clairement, il faut la réimplémenter ou sinon seule la "moitié" des coefficients sera prise en compte. On peut reprendre le travail effectué pour les classes mères, avec une légère modification. Rappelons que l'opération consistait à poser y = 0 puis à parcourir les éléments non nuls de la matrice en faisant à chaque fois y(i)+=a(i,j)*x(j). Il suffit donc de rajouter l'opération y(j)+=a(i,j)*x(i) lorsque i est différent de j.

Pour conclure cette section, notons que les modifications effectuées n'augmente pas la complexité en nombre d'opérations élémentaires. En effet, même si ce nombre augmente pour chaque élément, les éléments stockés sont deux fois moins nombreux : le nombre total reste le même.

3.2.5. Surcharge d'opérateurs pour employer la notation A(i,j)

Nous désirons maintenant "surcharger" l'opérateur () de façon ce que l'écriture A(i,j) désigne le coefficient de coordonnées (i,j) de la matrice A. La principale difficulté est que cette notation peut à la fois être utilisée pour :

• récupérer le coefficient, comme dans les expressions "A(3,9) < 2" ou "det = A(1,1)*A(2,2) - A(2,1)*A(1,2)"
• modifier le coefficient, comme dans les expressions "A(4,5) = 6" ou "A(1,2) += 3"

Nous remarquons alors qu'il faut à la fois surcharger l'opérateur () mais aussi l'opérateur = et éventuellement d'autres +=, -=, *=, /=, ++...

Une manière de réaliser cela est d'utiliser une classe Proxy représentant le coefficient $a_{i,j}$. Elle contient donc simplement une référence vers la matrice A et les coordonnées i, j. On commence alors par surcharger l'opérateur () à deux arguments des classes matrices pour qu'il retourne une instance de Proxy avec les arguments adéquats.

L'étape suivante est de surcharger des opérateurs sur la classe Proxy pour réaliser la lecture ou l'écriture. On surcharge l'opérateur () à zéro argument pour qu'il retourne A.getValue(i,j). Ainsi, lorsque l'on utilisera la notation A(i,j) dans une expression, on obtiendra le coefficient voulu. Pour permettre la modification, on surcharge l'opérateur =, de sorte qu'il effectue l'opération A.setValue(i, j). Nous pouvons généraliser cela aux autres opérateurs d'écriture et choisissons d'autoriser +=, -=, *= et /=. Il suffit de modifier légèrement setValue en lui ajoutant un paramètre optionnel indiquant la manière dont le coefficient doit être mis à jour (remplacement par value, addition de value...).

Remarquons qu'un objet Proxy est construit et détruit à chaque utilisation de l'écriture A(i,j). Cet objet est de taille négligeable : il ne contient qu'un pointeur (référence vers A) et deux entiers (les coordonnées i et j). A priori, cela ne devrait donc pas modifier beaucoup l'opération de remplissage. Vérifions-le avec la même matrice que dans la section 3.2.3. Pour des petites valeurs de la matrice le changement n'est pas perceptible. On a donc tracé le temps de remplissage en millisecondes pour des valeurs allant de 50000 à 3000000 (figure 6). On remarque que les droites sont pratiquement superposables, ce qui confirme que le coût supplémentaire est négligeable. Pour la plus grande matrice testée, il vaut 1 à 2% du coût initial.

Notons pour terminer qu'il aurait aussi été possible d'autoriser A[i][j] pour se rapprocher de la syntaxe du C++. Mais il aurait alors fallu passer par deux proxys intermédaires pour lire les coordonnées, ce qui aurait alourdit inutilement le code.

3.2.6. Importation et exportation de matrices

Les formats Harwell-Boeing (CSC) et Matrix-Market (COO) sont deux formats très utilisés pour stocker des matrices creuses sous forme de fichiers textes. Il existe des collections de matrices creuses (telles que [16] et [19]) intervenant dans divers domaines et applications. Permettre l'importation de ces matrices employées dans des problèmes concrets sera un bon moyen pour tester nos interfaces de solveurs. De plus, une option d'import/export est une fonctionnalité importante pour que les utilisateurs puissent échanger leurs matrices et comparer leurs résultats.

Par conséquent, j'ai ajouté aux classes de matrices macreuses correspondantes des méthodes offrant une telle possibilité. Pour cela j'ai utilisé deux librairies existantes :

• Matrix Market I/O library for ANSI C [1]
• C Routines for Harwell-Boeing File I/O [2]

Ces librairies ont été légèrement modifiées pour corriger quelques avertissements de compilation et permettre l'utilisation d'un paramètre générique pour les coefficients de matrices. J'ai aussi corrigé des fuites mémoires détectées dans la librairie HBIO et reportées par Valgrind.

Toutes les fonctionnalités mises à disposition ne sont pas utilisées dans les classes de matrices creuses. Ainsi, l'implémentation actuelle ne prend en compte ni les coefficients pattern/complexes ni les matrices denses.

La signature des méthodes est très simple. Elles prennent pour argument le nom du fichier à lire/écrire et la valeur de retour indique le succès ou non de l'opération :

• read_MatrixMarket et write_MatrixMarket. Elles sont disponibles dans les classes COO et COOsym. Elles acceptent uniquement les matrices creuses à coefficients réels ou entiers. Naturellement, COO traitent des matrices "générales" et COOsym des matrices "symétriques".
• read_HarwellBoeing et write_HarwellBoeing. Elles sont disponibles dans les classes CSC et CSRsym. Comme précédemment, CSC gère les matrices "RUA" (real unsymmetric assembled) tandis que CSRsym prend en compte des matrices "RSA" (real symmetric assembled).

3.2.7. Implémentation de sparse

Nous nous intéressons maintenant à l'implémentation d'une classe sparse<T> pouvant changer de type de stockage. Pour cela on considère quatre types sparse_COO, sparse_CSR, sparse_CSC et sparse_UNDEFINED ainsi que la propriété être symétrique ou non. L'utilisateur peut alors sélectionner le type grâce aux méthodes setSparseType(sparse_type) et setIsSymmetric(bool). En interne, une instance des classes COO, CSR, CSC, CSRsym ou COOsym est stockée sauf pour sparse_UNDEFINED.

Rappelons que dans le diagramme UML initial, sparse<T> implémente sparse_interface<T> et est composée de trois objets coo_, csr_, csc_. Ces objets sont alloués dynamiquement et libérés chaque fois que l'utilisateur change de format. Grâce au polymorphisme d'héritage, il est tout à fait possible d'utiliser les pointeurs coo_ et csr_ pour allouer des instances des classes COOsym et CSRsym. Pour que sparse implémente l'interface sparse_interface<T>, une méthode simple consiste à ce que chaque opération appliquée à l'objet sparse<T> utilise l'opération correspondante dans la classe actuellement instanciée. Il faut cependant prendre garde que les méthodes des classes parentes (setValue, getValue...) aient été déclarées virtual pour que la version des classes symétriques soient effectivement appelées. Lorsque la classe sparse<T> appelle les opérations friend, (A+B, k*A...), il est par contre nécessaire d'indiquer explicitement un cast statique. La plupart des opérations mettent en jeu une seule matrice et il est naturel de conserver le type de stockage. Le cas de l'addition/soustraction est particulier car les structures activées dans A et B ne sont pas forcément les mêmes. Nous partons du principe que lorsque que l'utilisateur effectue l'opération A += B, A conserve les mêmes structures actives. On traduit alors AplusB = A + B par AplusB = A; AplusB += B de sorte que dans une somme, le résultat prend les structures actives de l'opérateur de gauche.

La classe sparse<T> possède quelques méthodes supplémentaires par rapport à l'interface sparse_interface<T>. Nous avons déjà cités les méthodes pour changer le type de stockage. D'autres permettent de vérifier le stockage actuel comme isCOO(), isCSR(), isCSC(), isSymmetric(), sparseType(). La méthode clear(), vide la matrice et la remet au format sparse_UNDEFINED. Lorsque la matrice est en sparse_UNDEFINED, il est possible de modifier sa taille à l'aide de resize(). Les fonctions coo(), csr() et csc() retourne une copie de la matrice en utilisant la classe spécifiée. De même si la matrice est symétrique il est possible d'utiliser coosym() et csrsym().

Il existe plusieurs constructeurs de sparse<T>. Le constructeur par défaut est une matrice non symétrique de taille nulle et de type sparse_UNDEFINED. Des constructeurs prenant en argument une classe COO, CSR, CSC, CSRsym ou CSCsym permettent la conversion en sparse<T>. Un constructeur de copie et l'opérateur = sont aussi disponibles comme pour les autres classes implémentant sparse_interface<T>. Le constructeur à trois arguments est donné ci-dessous. Les valeurs par défaut correspondent au type CSR non symétrique. C'est le plus utilisé par les solveurs et il donne les meilleurs performances d'après la section 3.2.3.

sparse (size_type size, bool isSymmetric=false, sparse_type type=sparse_CSR)

Enfin, la classe sparse<T> dispose de méthodes pour lire et écrire une matrice. Rappellons que les types reconnus sont les suivants :

• MM_unsymmetric: MatrixMarket, format COO.
• MM_symmetric: MatrixMarket, format COO Symétrique.
• HB_unsymmetric: HarwellBoeing, format CSC.
• HB_symmetric: HarwellBoeing, format CSR Symétrique.

La signature des méthodes est données ci-dessous. Elles prennent pour premier argument le nom du fichier et retournent le succès ou non de l'opération. La fonction d'écriture ne fonctionne qu'avec les quatre types cités ci-dessus et le format est déterminé automatiquement. La fonction de lecture détruit le format actuel. Elle utilise le second argument pour déterminer le type du fichier. AUTODETECT_TYPE signifie qu'elle doit les essayer tous. Comme cela peut conduire à l'affichage d'erreurs lors des différentes tentatives, un troisième paramètre permet de les désactiver.

bool read (const char *filename, file_type type=AUTODETECT_TYPE, bool verbose=true)
bool write (const char *filename)

Pour terminer, notons que la classe possède des pointeurs privés pointer, subpointer et values sur les structures internes de la classe actuellement instanciée. Cela permettra aux classes de solveurs qui seront friend de la classe sparse<T> d'accèder aux données de la matrice.

3.2.8. Diagramme de classe révisé

Le travail effectué dans les sections précédentes nous permet d'apporter des modifications dans notre hiérachie de classes. Le diagramme révisé (figure 13) est disponible en Annexe D.

La principale différence est l'ajout de la classe sparse_common pour regrouper les implémentations de COO, CSR et CSC. Notons que les tableaux s'appellent maintenant pointer (auparavant rows dans COO), subpointer (columns dans COO/CSR, rows dans CSC) et values.

La classe Proxy, introduite dans la section précédente, est quant à elle utilisée par les classes implémentant l'interface sparse_interface<T>. Une telle classe est représentée par un paramètre M. En outre, toutes les classes utilisent le paramètre T pour les scalaires.

Deux types d'entiers ont été ajoutés. coord représente les coordonnées dans la matrice ou indices dans un tableau. size_type représente une taille ou un nombre d'éléments. Ils étaient initialement implémentées comme des unsigned int. Mais ils ont été ensuite redéfinis comme des entiers signés car tous les solveurs utilisent des int.

3.2.9. Iterateurs sur les matrices

Cette option a été ajouté après étude des besoins d'H2OLAB (cf section 4.9) à la fin de mon stage. Par conséquent, je n'ai eu le temps que d'implémenter les itérateurs pour la classe sparse<T>. De plus les opérations d'édition sont autorisées même quand la matrice est de type const.

L'idée est de fournir un mécanisme de parcourt similaire à celui que l'on trouve pour des objets de la librairie standard tels que vector<T>. On utilise pour cela une classe sparse<T>::iterator dont l'instance va servir de curseur sur la matrice. Une telle instance est crée par un appel aux membres iterator_begin() ou iterator_end() plaçant respectivement au début et à la fin de la matrice. Notons que pour CSR (respectivement CSC) il est possible d'indiquer comme argument le numéro de ligne (respectivement colonne) sur lequel on veut placer l'itérateur. La classe sparse<T>::iterator utilise les même techniques que sparse<T> pour améliorer la syntaxe : surchage d'opérateurs et utilisation d'un proxy ij_value. Ainsi la boucle ci-dessous multiplie tous les éléments diagonaux par 10 et affiche les valeurs des éléments extra-diagonaux :

for (sparse<T>::iterator it = A.iterator_begin(); it != A.iterator_end(); it++)
  {
     if (it->isDiagonalElement())
       *it *= 10;
     else
       cout << "A(" << it->i() << "," << it->j() << ") = << " *it << endl;
  }

Les opérations disponibles dans la classe sparse<T>::iterator sont :

• *it qui retourne la valeur $a_{i,j}$ pointée par l'itérateur.
• it->i() et it->j() qui retournent les coordonnées de l'itérateur dans la matrice.
• it->isDiagonalElement() qui indique si l'itérateur est sur un élement diagonal.
• *it = v, *it += v, *it -= v, *it *= v, *it /= v qui modifie l'élément $a_{i,j}$ de la matrice.

4. Hiérarchie de classes pour les interfaces de solveur

4.1. Implémentation de quelques solveurs séquentiels

Avant d'aborder l'interface de solveurs parallèles, il m'a été proposé d'implémenter quelques solveurs séquentiels. Cela m'a permis de tester la classe sparse_interface<T> et de dégager les premières idées de paramètres et statistiques communes aux solveurs.

Considérons tout d'abord les points communs des solveurs. Ceux-ci seront si possible regroupés dans une classe solver. Chacun cherche à résoudre un système $Ax=b$ , où A est une matrice de type sparse_interface<T> et b est un vector<T> de même taille. Le constructeur de l'objet prend donc en compte ces deux paramètres. Ceux-ci sont en fait des références const de façon à éviter d'avoir à recopier ces objets qui peuvent prendre beaucoup de place en mémoire.

Une fonction solve() lance l'exécution de la recherche de solution et stocke le résultat x. Ce dernier peut alors être récupéré par une fonction membre x(). Ajoutons deux fonctions membres intéressantes : residual_norm() calculant la norme résiduelle $\frac{\parallel \mathrm{Ax}-b\parallel }{\parallel b\parallel }$ (précision de la solution) et time() donnant le temps d'exécution qui a été nécessaire pour déterminer cette solution (efficacité en temps).

Ensuite, j'ai choisi quelques méthodes classiques : Jacobi, Gauss-Seidel, SOR et le Gradient Conjugué. Les algorithmes pour ces méthodes étant bien connus, la principale difficulté est ici de les adapter pour le cas des matrices creuses. Rappelons qu'il s'agit de méthodes itératives basées sur le même principe : on part d'un vecteur initial $x^0$ et à l'étape k+1, on déduit un vecteur $x^{k+1}$ à partir de $A,b,x^k$ . On espère alors que la suite ${\left(x^k\right)}_{k\in \mathbb{N}}$ converge vers une solution x. On considère deux critères d'arrêts :

• la norme résiduelle est inférieur à une certaine valeur de tolérance ε
• on a atteint le nombre d'itérations maximum ${\mathrm{iter}}_{max}$

Cela nous conduit à trois méthodes qui seront communes à une certaine classe IterativeSolver : set_tolerance(), set_max_iterations() et iterations(). Les deux premières ont une signification claire, la dernière permet de connaitre le nombre d'itérations effectuées. Il peut aussi être important de fixer la valeur initiale du sytstème. Par défaut, $x^0=0$ mais la méthode set_initial_guess() permet de changer cette valeur. Si aucun paramètre ne lui est passé, il met x à 0. Sinon, il met x à la valeur donnée.

Pour les méthodes de Jacobi, Gauss-Seidel et SOR, les coordonnées du vecteur $x^{k+1}$ sont données par :

• ${x_i}^{k+1}=\frac{1}{a_{i,i}}(b_i-\sum _{j\ne i}{a}_{i,j}{x_j}^k)$ pour Jacobi
• ${x_i}^{k+1}=\frac{1}{a_{i,i}}(b_i-\sum _{j<i}{a}_{i,j}{x_j}^{k+1}-\sum _{j>i}{a}_{i,j}{x_j}^k)$ pour Gauss-Seidel
• ${x_i}^{k+1}=\left(1-\omega \right){x_i}^k+\frac{\omega }{a_{i,i}}(b_i-\sum _{j<i}{a}_{i,j}{x_j}^{k+1}-\sum _{j>i}{a}_{i,j}{x_j}^k)$ pour SOR de paramètre $\omega$

Par conséquent, pour chaque i, il faut parcourir toute la ligne i de A et déterminer le coefficient $a_{i,i}$ . Rappelons cette remarque pour Gauss-Seidel et SOR : le coefficient ${x_i}^{k+1}$ n'utilise pas les valeurs de ${x_j}^k$ pour j < i. Autrement dit dans le tableau où x est stocké, on peut directement écraser ${x_i}^k$ avec la valeur de ${x_i}^{k+1}$ si on parcourt les lignes i dans l'ordre. Par contre pour Jacobi, le passage à une variable temporaire x_next est nécessaire.

Pour obtenir les coefficients $a_{i,j}$ on n'utilise pas la méthode getValue, qui demanderait à chaque fois de parcourir les premiers éléments d'une ligne. Il est donc nécessaire de mettre les classes Jacobi, GaussSeidel et SOR en friend de sparse<T> pour pouvoir accéder directement aux structures de la matrice creuse. Les considérations précédentes conduisent à choisir le format CSR pour la matrice A : on utilise une variable i parcourant les lignes de la matrice A et une variable k parcourant la ligne i, i.e. allant de A.pointer[i] à A.pointer[i+1] - 1. Il faut simplement déterminer à chaque étape i l'(unique) indice k_diag vérifiant A.column[k_diag] == i. Le coefficient $a_{i,i}$ utilisée dans les formules de récurrence est alors A.values[k_diag]. Si l'indice k_diag n'est pas trouvé, alors la matrice A comporte un élément nul sur la diagonal et la précondition d'utilisation des méthodes n'est pas satisfaite.

La méthode du gradient conjugué est une méthode vectorielle et par conséquent ne nécessite pas la connaissance des structures internes des matrices creuses. C'est pourquoi la matrice sparse<T> utilisée peut avoir n'importe quel format. On implémente ensuite dans le fichier vector_utils.h quelques opérations élémentaires :

• somme et différence de vecteurs : u+v, u-v
• norme euclidienne d'un vecteur : norm(u)
• produit scalaire de deux vecteurs : u*v
• multiplication d'un vecteur par un scalaire : k*u

Ces opérations sont pour la plupart surchargées de façon à rendre plus lisible l'algorithme du Gradient Conjugué. Combinées avec la multiplication matrice vecteur, elles sont suffisantes pour implémenter cet algorithme. Notons aussi qu'elles peuvent être utiles ailleurs, par exemple pour obtenir la norme résiduelle : norm(A*x - b) / norm(b).

4.2. Analyse générale des solveurs à interfacer

Les solveurs pour lesquels une interface m'a été demandée sont Hypre [11], MUMPS [15], SuperLU disttribué [18], SuperLU [17] et UMFPACK [20]. Les trois premiers sont des solveurs parallèles tandis que les deux derniers sont des solveurs séquentiels. Il est alors intéressant d'introduire deux classes ParallelSolver et SequentialSolver héritant de la classe Solver, pour implémenter les caractéristiques propres (figure 7).

Rappellons aussi que nous avons introduit une classe IterativeSolver dans la section précédente. Nous tirons profit de l'héritage multiple du C++ pour que Jacobi, GaussSeidel, SOR et ConjugateGradient héritent de cette classe en plus de SequentialSolver. De même, Hypre regroupant des méthodes itératives, nous le faisons aussi hériter de cette classe. Finalement, nous obtenons l'héritage représenté sur la figure 8.

Remarquons toutefois que IterativeSolver doit pouvoir accéder à la variable x pour pouvoir modifier la valeur initiale dans set_initial_guess(). Par conséquent, elle doit aussi hériter de la classe Solver. Nous sommes alors conduit au "problème du losange". Par exemple Jacobi hérite à la fois de SequentialSolver et IterativeSolver qui héritent eux-même tous les deux de la classe Solver. Par défaut, le langage C++ va dupliquer les variables membres de Solver et il y aura une ambiguïté lorsqu'on voudra les utiliser dans la classe Jacobi. Le même problème se pose pour GaussSeidel, SOR, ConjugateGradient et Hypre. Pour remédier à cela, on déclare les classes SequentialSolver, ParallelSolver et IterativeSolver avec un héritage virtuel. Ainsi, les objets de Solver ne seront pas dupliqués. La hiérarchie générale des classes de solveurs est donnée en Annexe D dans le diagramme de la figure 14.

4.3. Listes des solveurs

Le tableau de la figure 9 résume les différents types de solveurs et les formats acceptés dans l'implémentation actuelle. L'utilisateur doit s'assurer qu'un format adéquat est donné au solveur.

4.4. La classe solveur parallèle

Cette classe regroupe les solveurs Hypre, SuperLU distribué et MUMPS. L'idée est de découper la matrices et les vecteurs en morceaux que l'on répartit sur plusieurs processeurs. Ceux-ci effectuent en parallèle les opérations nécessaires à la résolution du système. Chacun fournit alors un morceau de la solution x. A la fin, on doit récupérer tous ces morceaux pour obtenir la solution du système.

MUMPS se distingue en demandant une saisie au format COO. La documentation ne donne pas de contraintes sur la manière donc la découpe de la matrice est réalisée et il est donc a priori possible de répartir un nombre (pratiquement) équitable d'éléments non nuls pour chaque processeur à l'instar de ce qui va être fait ci-dessous. En fait, MUMPS autorise deux types de saisie : soit on envoit chaque morceau sur un processeur, soit on donne l'intégralité de la matrice sur un processeur hôte. Comme la matrice sparse<T> a été construite sur un processeur, on utilise la saisie centralisée et on laisse à MUMPS le soin de répartir les données sur les processeurs.

Les autres solveurs utilise un format d'entrée CSR ou CSC. Les variables à fournir varient beaucoup selon les solveurs mais suit le schéma général représenté sur la figure 10. Chaque processeur $p_k$ reçoit un groupement de n_loc lignes/colonnes contiguës à partir d'une certaine valeur de départ i_lower. L'ensemble forme une partition de la matrice. Le découpage obtenu coïncide avec celui du second membre b et de la solution x.

Une façon simple de répartir les données est d'essayer de donner un même n_loc à tous les processeur. Citons deux possibilités :

• On calcul la valeur $\lfloor \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.\rfloor$ qui est le nombre de lignes/colonnes reçues par chaque processeur, sauf le dernier qui prend ce qui reste. On a $\lfloor \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.\rfloor \le \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.$ donc $n-(\mathrm{nbproc}-1)\lfloor \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.\rfloor \ge \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.\ge \lfloor \raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$\mathrm{nbproc}$}\right.\rfloor$ . Ainsi le dernier processeur reçoit toujours plus d'éléments que les autres.
• On réalise la division euclidienne de $n=\mathrm{nbproc}q+r$ , avec $r<\mathrm{nbproc}$ . Chaque processeur reçoit alors q éléments. On répartit les r éléments restants parmi r processeur.

Une approche alternative est d'essayer de répartir un nombre équitable d'éléments nnz dans chaque processeur. A cause de la contrainte sur les lignes, ce n'est bien sûr pas possible mais on essaye de s'en approcher. Pour cela, on commence par calculer un nnz_moyen = arrondi(nnz / nbproc). L'idée est de parcourir les lignes/colonnes et de regrouper des lignes dans un même processeur, de façon à ce que leur nombre d'éléments non nuls soient le plus proche possible de nnz_moyen. Si une ligne comporte plus de nnz_moyen éléments, alors on a pas le choix on doit mettre ces éléments entièrement sur un processeur. Sinon on regroupe les lignes de façon à ne pas dépasser le nombre nnz_moyen. On préfère donner une limite supérieure car on espère que cela va s'équilibrer avec le cas des processeur qui reçoivent une seule ligne contenant beaucoup plus d'éléments non nuls que la moyenne.

Pour éviter que certains processeur ne reçoivent aucune ligne/colonne, on s'assure aussi que le nombre de lignes/colonnes restants est suffisant. Cela entraine une dissymétrie dans l'algorithme car il se peut alors qu'on ne cherche plus à assurer l'équilibre du nombre d'élément pour les derniers processeur. Toutefois il s'agit ainsi d'un choix simple évitant d'avoir recourt à un retour en arrière. Notons aussi que s'il y a plus de processeur que de lignes/colonnes dans la matrice, certains processeur ne recevront aucune ligne. Heureusement cela ne se produit presque jamais en pratique.

Notons que dans le cas parfaitement symétrique où toutes les lignes/colonnes ont le même nombre d'éléments non nuls et où nbproc divise le nombre de ligne/colonne alors les trois algorithmes précédents donnent la même répartition. Si la complexité des algorithmes de résolution dépendent du nombre d'éléments non nuls alors essayer de répartir ce nombre d'éléments plutôt que le nombre de lignes semble une meilleur idée. L'algorithme Scatter est donné en Annexe B

Nous n'avons considéré ici que la répartition des lignes/colonnes de la matrice. Dans le cas d'une matrice carrée, cette répartition coïcide avec celle des lignes de la solution et du second membre. L'algorithme précédent est exécuté sur le processeur hôte. Celui-ci conserve les parties qu'il va traiter et envoie les autres parties au processeur concerné. Les autres processeur se contentent de recevoir les données envoyées par le processeur hôte. A la fin de la résolution, une fonction gather réalise l'opération inverse : chaque processeur envoie la partie de la solution au processeur hôte. La classe ParallelSolver comporte enfin quelques fonctions supplémentaires pour initialiser les paramètres MPI ou encore mesurer le temps de calcul en prenant le maximum du temps d'exécution de chaque processeur.

Pour conclure cette section, remarquons que dans tous les solveurs parallèles, on fait l'hypothèse que les données d'entrée/sortie sont stockées sur le processeur de rang zéro. Par conséquent, l'utilisateur doit bien prendre garde de construire la matrice/vecteurs et récupérer les résultats sur ce processeur. Construire la matrice sur les autres processeurs est inutile et entraine une perte de temps, surtout si la matrice est très grosse. La validité des résultats est quant à elle assurée uniquement pour le processeur zéro, les valeurs sur d'autres processeurs ne sont pas définies.

4.5. Spécificité de chaque solveur

En plus des fonctions membres qui ont été présentées plus haut, les solveurs possèdent leur propres fonctions de paramètrage (de la forme set_*) et de récupération de statistiques (de la forme get_*). Les solveurs à interfacer possèdent un nombre très important d'options et de résultats. Implémenter l'intégralité des possibilités offertes par les solveurs aurait été fastidieux et non forcément utile pour l'utilisateur lambda. De plus, cela risquerait d'alourdir l'interface. Ainsi il a été fixé comme objectif minimum d'implémenter les fonctionnalités déjà utilisées par l'équipe SAGE, sachant que l'interface pourrait être augmentée si nécessaire.

Chaque solveur comporte un format de prédilection qui devrait être utilisé lors de l'appel au constructeur (cf section 4.3). Ainsi Hypre, SuperLU et SuperLU_dist utilisent le format CSR, MUMPS le format COO et UMFPACK le format CSC.

UMFPACK comporte trois étapes : factorisation symbolique, factorisation numérique et résolution. Comme on cherche à fournir une interface simple, la fonction solve se contente de réaliser les trois étapes en une fois. Le temps mis pour chacune des étapes est toutefois mesuré et peut être récupéré par l'utilisateur. UMFPACK autorise un petit nombre de paramètres de contrôle qui sont contenus dans un tableau. Dans l'interface proposée, l'utilisateur n'a pas besoin de remplir ce tableau à la main, mais doit utiliser des fonctions membres set_* avec un nom plus explicite.

MUMPS possèdent lui aussi un certain nombre de paramètres fournis à travers des tableaux ICNTL (pour les entiers) et CNTL (pour les scalaires). Des fonctions avec des noms explicites permettent une configuration plus compréhensible. Toutefois, des fonctions set_icntl(int index, int value) et set_cntl(int index, double value) permettent aux utilisateurs avancés de modifier les paramètres non accessibles avec l'interface standard. Comme pour UMFPACK, il est possible de récupérer le temps mis pour chaque étape (analyse, factorisation et résolution) ainsi que quelques statistiques supplémentaires.

Les solveurs SuperLU et SuperLU_dist possèdent des paramètres et statistiques communes. Une classe superlu_common regroupe les fonctions membres permettant de gérer ces données. Les classes SuperLU et SuperLU_dist héritent alors de cette classe, ce qui leur permet de partager une partie de leur implémentation. Notons que SuperLU accepte le format CSR ou CSC alors SuperLU_dist ne prend en compte que le CSR. Pour rester cohérent, j'ai retenu le format CSR pour les deux solveurs. Enfin, SuperLU_dist possèdent une particularité : il répartit les processeur selon une grille. La fonction set_grid_size (int nprow, int npcol) permet à l'utilisateur de spécifier une répartition des processeur selon un rectangle de taille nprow et npcol. Par défaut, le choix qui est fait est une grille carrée de dimension $\lfloor \sqrt{\mathrm{nbproc}}\rfloor$. L'utilisateur doit faire en sorte que nbproc = npcol * nprow : la résolution échoue s'il n'y a pas assez de processeur (nbproc < npcol * nprow) et les processeur en trop sont gaspillés s'il n'y a pas assez de place sur la grille (nbproc > npcol * nprow).

Enfin Hypre est constitué d'un ensemble de solveurs et préconditionneurs. La fonction set_solver (hypre_solver solver, hypre_solver preconditioner=NONE) permet à l'utilisateur de choisir un couple solveur/preconditionneur. Les associations retenues sont les suivantes :

• BiCGSTAB en solveur, avec un préconditionneur facultatif parmi BoomerAMG, ParaSails ou Euclid.
• BoomerAMG utilisé comme solveur.
• GMRES en solveur, avec un préconditionneur facultatif parmi BoomerAMG, ParaSails ou Euclid.
• PCG en solveur, avec un préconditionneur facultatif parmi BoomerAMG, ParaSails ou Euclid.

Notons que les fonctions de paramètrage/statistiques héritées de iterative_solver s'applique au solveur. Il est possible de fixer la tolérance ou le nombre maximum d'itérations du preconditionneur avec des fonctions set_preconditioner_*. Par défaut, elles sont mises respectivement à 0 et 1, de façon à ce que le préconditionneur s'applique exactement une fois. Pour distinguer les options propres à chaque solveur/préconditionneur, les fonctions correspondantes s'écrivent sous la forme set_NomDuSolveurOuPreconditionneur_*.

4.6. Préconditionnement et mise à l'échelle des matrices

Il est souvent utile de multiplier la matrice (et le vecteur) du système par d'autres matrices pour accélérer la résolution du système. Des opérations particulièrement simples sont les mises à l'échelle (multiplication par une matrice diagonale) ou les permutations des lignes/colonnes (multiplication à gauche ou à droite par une matrice de permutation). Je n'ai pas travaillé sur l'intégration de librairie comme HSL [10] permettant des mises à l'échelle.

Il existe divers algorithmes pour déterminer des permutations efficaces. Certains sont implémentés dans les solveurs à interfacer d'autres dans des librairies indépendantes. Il est difficile de trouver une interface commune. Idéalement, il faudrait permettre de préconditionner les matrices à l'extérieur des interfaces de solveurs, mais cela peut poser des problèmes de gestion du second membre et de la solution. Actuellement, je me suis contenté d'implémenter les méthodes à l'intérieur des fonctions solve de chaque solveur.

• MUMPS supporte des permutations à appeler avec la fonction set_ordering. Il est aussi possible d'utiliser mumps_permuting_scaling ou mumps_scaling_strategy, même s'il vaut mieux dans le doute laisser les valeurs par défaut qui demande à MUMPS de rechercher la meilleur stratégie.
• SuperLU et SuperLU_dist possèdent une fonction set_permutation_methods(superlu_permute_row permute_row, superlu_permute_col permute_col).

4.7. Expérimentation des interfaces de solveurs

J'ai compilé et testé mon programme sur ma machine locale ainsi que sur les machines de la grille. Les solveurs parallèles ont été testés avec 2, 4, 8 ou 16 processeurs. Les principales commandes à réaliser sont les suivantes (N et T sont des paramètres à remplacer) :

login@localhost% ssh login@paramount                        # connection sur la machine frontale
login@paramount% oarsub -I -l nodes=N,walltime=T            # réservation de N noeuds pendant une durée T
login@paradent-x% mpdboot -n N -f $OAR_NODEFILE --rsh=oarsh # Lancement du démon MPI sur les N machines réservées
login@paradent-x% mpiexec -np N ./test_parallel_solver.o    # exécution d'un programme de test sur les N machines

Pour vérifier qu'un solveur fournit une bonne solution pour une matrice A donnée, on construit un vecteur x (par exemple $x={\left(1,2,...,n\right)}^t$) et on génère le second membre $b=Ax$. Il faut alors s'assurer que la solution calculée par le solveur correspond bien au x initialement construit. Pour les solveurs directs, j'ai vérifié que la norme résiduelle (fonction residual_norm()) était très faible. Pour les solveurs itératifs, elle doit être en dessous du seuil de tolérance si le nombre d'itérations maximum n'a pas été atteint. Certaines matrices nécessitent de bien paramètrer le solveur pour obtenir une solution convenable. Des cas de tests m'ont été fournis par Désiré Nuentsa. Diverses matrices que j'ai testées ainsi que leurs caractéristiques sont données dans le tableau de la figure 11.

Nom	Taille (n)	Nombre d'éléments non nuls (nnz)	Propriété	Solveurs utilisés	Aperçu de la structure
mbeause	496	41 603	non symétrique	Hypre (GMRES)
sherman2	1 080	23 094	non symétrique	Hypre (GMRES+Euclid)
sherman4	1 104	3 786	non symétrique	Hypre (BiCGSTAB), MUMPS, SuperLU_dist, UMFPACK, SuperLU
bcsstk14	1 806	32 630	symétrique	Hypre (BiCGSTAB), MUMPS
orsreg_1	2 205	14 133	symétrique	Hypre (BiCGSTAB), MUMPS, SuperLU_dist, UMFPACK, SuperLU
mark3jac020sc	9 129	52 883	non symétrique	MUMPS et SuperLU_dist
bbmat	38 744	1 771 722	non symétrique	Hypre (GMRES+Parasails)
s3dkq4m2	90 449	2 455 670	symétrique	Hypre (BiCGSTAB), MUMPS, SuperLU_dist, UMFPACK
nlpkkt160	8 345 600	225 422 112	symétrique	Hypre (GMRES)

figure 11

J'ai aussi comparé le temps de résolution et la précision des deux méthodes de découpage vues dans la section 4.4 : nombre de lignes/colonnes ou nombre d'éléments non nuls équitablement répartis. Les tests ont été menés sur certaine des matrices précédentes. A cause de la symétrie des problèmes, ces matrices ont souvent un nnz par ligne régulier. Ainsi, j'ai aussi utilisé une matrice issue de l'économie "mbeause" avec un nnz par ligne très variable ainsi qu'une matrice artificielle en forme de triangle (nnz par ligne croissant). Aucun des tests effectués n'a montré une amélioration significative des performances lors de la répartion équitable du nombre d'éléments non nuls. Il semble donc préférable de conserver une répartition équitable des lignes car l'algorithme est plus simple.

4.8. Limites de l'implémentation

Citons quelques limites de l'implémentation actuelle :

• La lecture Harwell-Boeing peut échouer pour certaines matrices, par exemple la "s3dkq4m2" du site de MatrixMarket [16]. En effet, la librairie utilisée [2] n'est pas capable de reconnaitre tous les descripteurs de format FORTRAN lorsque ceux-ci sont trop complexes.
• Le paramètre générique T des matrices permet de manipuler des float ou double pour que l'utilisateur puisse jouer sur la précision ou la place en mémoire. Néanmoins les interfaces des solveurs utilisent des double et c'est donc celle qui est effectivement utilisées. UMFPACK autorise aussi des nombres complexes, mais aucune interface n'a été mise en place pour cette option. Par conséquent l'implémentation actuelle ne permet pratiquement que d'utiliser des double . Il serait nécessaire d'apporter quelques modifications si on souhaite généraliser l'utilisation à d'autres type T .
• Certains solveurs autorisent le passage de plusieurs seconds membres b. Cette possibilité n'a pas été retenue dans l'interface car elle n'était pas utilisée par l'équipe SAGE.
• Le préconditionnement des matrices devrait être amélioré et complété (cf section 4.6).
• La compilation d'un programme utilisant Hypre et SuperLU_dist entraine des conflits de nom. Il faut alors autoriser les redéfinitions.
• Il n'est pas possible de compiler un programme utilisant simultanément SuperLU et SuperLU_dist.

4.9. Comparaison avec les implémentations existantes

L'implémentation actuelle reprend les différentes paramètres et statistiques de l'interface C réalisée par Désiré Nuensta. La structure générale du programme a été grandement améliorée grâce à la conception sous forme de classes. Cela a notamment permis de simplifier le code en partageant les points communs des solveurs pour qu'ils puissent être manipulés de manière uniforme. Pour les méthodes de préconditionnements (cf section 4.6), toutes celles du code original n'ont pas été implémentées/testées. Certains formats d'entrées/sorties ou d'interface de solveurs encore en développement n'ont pas été intégrées. Le travail effectué donne toutefois une base solide pour des extensions futures.

A la fin de mon stage, j'ai aussi analysé le code d'H2OLAB et préparé l'intégration de mon module. Une première remarque est que H2OLAB utilise une structure de maps [23] pour représenter ses matrices. Les clés sont les coordonnées (i, j) et les valeurs correspondantes les $a_{i,j}$. L'étude du code montre que le format CSR semble bien adapté pour remplacer ces maps, puisque les algorithmes lisent "en ligne". Notons toutefois quelques régressions :

• H2OLAB utilise le solveur SMG de Hypre, qui ne peut être utilisé qu'avec l'interface "Structured-Grid System". Mais celle-ci ne semble pas pouvoir être construite avec une matrice creuse.
• H2OLAB peut utiliser UMFPACK avec des nombres complexes. Cela n'est actuellement pas possible avec la classe sparse<T>.
• Il est nécessaire de préciser le début et la fin de l'édition de la matrice pour la classe sparse<T>.

Par contre le format sparse<T> en CSR pourrait théoriquement apporter des améliorations par rapport aux les maps :

• L'implémentation des maps nécessitent une structure d'arbre binaire pour ranger les valeurs. Cela utilise plus de place en mémoire que les tableaux au format CSR.
• Contrairement aux maps, les tableaux de CSR sont stockées de manière contigüe. Comme les caches des processeurs tiennent compte de la localité spatiale, il se pourrait qu'appliquer des opérations sur les éléments de la matrice soit plus rapide avec CSR.
• La récupération d'un élément se fait en O(log(nnz)) pour les maps au lieu de O(log(nnz par ligne)) pour CSR.
• Pour le remplissage, chaque insertion se fait en O(log(nnz))+"temps d'allocation d'un élément" pour une map. Pour CSR, les tableaux ne sont redimensionnés qu'au début et à la fin de l'édition. Si le slot_size a bien été choisi, les éléments sont directements insérés en O(log(nnz par ligne)). Un choix de slot_size semble être possible pour H2OLAB du fait de la régularité des problèmes physiques traités.
• Pour l'envoi de la matrice au solveur, CSR peut directement être envoyé au solveur. Si cela est nécessaire la matrice est convertie au format COO/CSC en O(nnz + size). A l'inverse, la map nécessite plusieurs remplissages de matrice : celui de la map et celui des solveurs utilisés.

L'étude du code source d'H2OLAB a aussi montré deux opérations manquantes pour sparse<T>: la suppression d'un élément et le parcourt de la matrice. La première opération consiste à enlever l'élément non nul à la position (i, j). Elle pourrait théoriquement être remplacée par un setValue(i, j, 0) mais alors l'élément serait simplement mis à zéro et non réellement retiré. Cela pourrait diminuer les performances puisque l'on stocke alors des éléments inutiles. Par conséquent, j'ai ajouté à l'interface sparse_interface<T> une opération remove(i, j) qui peut être appelée durant l'édition d'une matrice. Elle retire l'élément si elle le trouve dans les tableaux internes ou dans la map toInsert.

La seconde opération, qui a été décrite dans la section 3.2.9, est le parcourt de la matrice. Très fréquemment, seule une ligne de la matrice est parcourue en excluant éventuellement l'élément diagonal. On ne peut alors plus utiliser l'itérateur sur une map. Typiquement le parcourt est réalisé dans H2OLAB de la façon suivante :

1. Énumération de toutes les coordonnées j des éléments extradiagonaux de la ligne i grâce à une fonction std::vector<int> one_equal_to(i). Cette fonction cherche la position d'un élément (i, i) et fait deux recherches avec des itérateurs en avant/arrière en partant de cette position, jusqu'à sortir de la ligne i. Les indices j des éléments trouvés sont stockés dans le vector<int> retourné par la fonction .
2. On parcourt le vector<int> des j et pour chacun on effectue un appel à une opération de la map accèdant à l'élément (i, j). Cet accès nécessitant à chaque fois de parcourir la map.

On constate que cette implémentation est très coûteuse. De plus l'étape 1 ne fonctionne que si l'élément (i, i) existe. On doit aussi à chaque fois réaliser un cas particulier pour cet élément diagonal selon que l'on veut ou non l'inclure dans l'étape 2. L'utilisation de sparse<T> au format CSR permet de simplifier considérablement ce parcourt en faisant simplement une boucle du type :

for (sparse<T>::iterator it = A.iterator_begin(i); it != A.iterator_end(i); it++)
  {
    if (it->isDiagonalElement())
      // opération sur *it
  }

Ici, les accès iterator_begin(i) et iterator_end(i) se font via pointer[i] et pointer[i+1] donc sont directs. Le curseur it permet aussi d'accèder directement à l'élément pointé. Enfin, une fonction isDiagonalElement() est disponible pour rendre le code plus lisible.

Enfin, la construction des matrices dans H2OLAB se fait en parallèle, c'est-à-dire qu'une map est construite sur chaque processeur. Ces maps comportent chacune une partie des éléments non nuls de la matrice. Cela est nécessaire pour rendre efficace le remplissage de la matrice et éviter qu'à un moment donné, l'intégralité de la matrice soit stockée sur un processeur ce qui peut poser des problèmes pour les grosses matrices. En ce qui concerne la classe sparse<T>, l'utilisateur peut actuellement construire les morceaux de la matrice en parallèle en utilisant des setValue distincts selon les processeurs.

A titre d'exemple, pour la construction d'une matrice diagonale :

A.begin() // executé par tous les processeurs
for (coord i = 0; i < SIZE; i++)
  {
    if (mpi_rank == (i%mpi_size))
      A.setValue(i, i, i); // exécuté par le processeur congru à i modulo mpi_size
  }
A.end() // executé par tous les processeurs

Néanmoins, ces morceaux doivent être répartis selon le découpage décrit dans la section 4.4. Pour résoudre ce problème, il est nécessaire de modifier les classes sparse_common et/ou parallel_solver. J'ai commencé à modifier la fonction end() de sparse_common pour que les morceaux soient redistribués à la fin de la construction. Il faut ensuite modifier sparse_common pour ne plus opérer de découpage lorsque la matrice a été construite en parallèle.

5. Conclusion

Le travail qui m'était demandé se décompose en deux parties, avec pour chacune une phase de conception et de développement :

• Réalisation d'une hiérarchie de classe pour manipuler différents formats de matrices carrées creuses.

  Les formats disponible sont COO et CSR/CSC. Il est aussi possible de manipuler des matrices symétriques ainsi qu'un objet sparse<T> de gestion de plusieurs formats. Les opérations de bases comme la somme, la transposée, les conversion ou le produit matrice-vecteur ont été implémentées de manière très efficace. Une technique spéciale a été utilisée pour le remplissage de matrice, de façon à éviter la perte de temps due à la lourdeur de l'opération "insertion d'un élément dans un tableau". J'ai réalisé une étude des performances de ces diverses opérations. J'ai aussi mis à disposition des fonctions d'écritures/lectures des formats de fichier MatrixMarket et HarwellBoeing.

• Réalisation d'interfaces vers des solveurs en utilisant le type matrice creuse précédent.

  Des solveurs séquentiels classiques Jacobi, Gauss-Seidel, SOR et Gradient Conjugué ont tout d'abord été implémentés. Ensuite, j'ai réalisé une hiérarchie de classes pour interfacer différents solveurs : UMFPACK, SuperLU, Hypre, SuperLU_dist, MUMPS. Cette hiérarchie a été fait de façon à ce que le maximum de partie commune aux différents solveurs partage la même implémentation. Diverses fonctions de paramètrage et de récupérations de résultats ont été mises en place pour les solveurs. Une attention particulière a été donnée aux solveurs parallèles, pour lesquels il est nécessaire de répartir les données sur plusieurs processeurs à l'aide du langage MPI. Des tests ont été effectués sur des grilles de calculs avec des matrices de différentes formes, de taille $n$ allant jusqu'à 8 millions et un nombre d'éléments non nuls nnz pouvant atteindre plusieurs millions. Les solveurs parallèles ont bien sur été testés avec plusieurs processeurs (2, 4, 8 ou 16).

Ce stage m'a donné la possibilité de mettre en pratique les connaissances de C++, de conception et de programmation orientée objet vue en seconde année à l'ENSIIE. J'ai pu notamment me servir des notions de patrons de classe, de polymorphisme d'héritage, d'héritage multiple, d'héritage virtuel... J'ai aussi pu utiliser les concepts d'analyse numérique de première année dans le cas de système creux de taille très élevée. J'ai pu entrevoir le fonctionnement de divers solveurs pour résoudre ce type de système. J'ai enfin découvert le langage MPI et les programmes exécutés en parallèle. J'ai appris à coder et débogguer ce type de programmes et testé leur exécution sur des grilles de calculs.

6. Bibliographie

1. ANSI C library for Matrix Market I/Ohttp://math.nist.gov/MatrixMarket/mmio-c.html
2. C Routines for Harwell-Boeing File I/O - http://math.nist.gov/~KRemington/harwell_io/harwell_io.html
3. Cmake, Cross Platform Make - http://www.cmake.org/
4. Doxygen - http://www.doxygen.org
5. GCC, the GNU Compiler Collection - http://gcc.gnu.org/
6. GDB: The GNU Project Debugger - http://www.gnu.org/software/gdb/
7. GNU emacs - http://www.gnu.org/software/emacs/
8. Grid'5000 - https://www.grid5000.fr/
9. H2OLAB - http://hydrolab.irisa.fr/
10. HSL - http://www.aspentech.com/hsl/
11. Hypre - http://acts.nersc.gov/hypre/
12. MPI-2 - http://www.mpi-forum.org/docs/mpi-20-html/mpi2-report.html
13. MPICH2 - http://www.mcs.anl.gov/research/projects/mpich2/
14. MTL4 manual, matrix Insertion - http://www.osl.iu.edu/research/mtl/mtl4/doc/matrix_insertion.html
15. MUMPS : a parallel sparse direct solver - http://mumps.enseeiht.fr/
16. Matrix Market - http://math.nist.gov/MatrixMarket/
17. SuperLU - http://crd.lbl.gov/~xiaoye/SuperLU/#superlu
18. SuperLU_dist - http://crd.lbl.gov/~xiaoye/SuperLU/#superlu_dist
19. The University of Florida Sparse Matrix Collection - http://www.cise.ufl.edu/research/sparse/matrices/
20. UMFPACK - http://www.cise.ufl.edu/research/sparse/umfpack
21. Umbrello UML Modeller - http://uml.sourceforge.net/
22. Valgrind - http://valgrind.org/
23. Wikipedia, article "map (C++ container)" - http://en.wikipedia.org/wiki/Map_(C++_container)

Annexes

A. Algorithmes de conversion et transposition

Fonction CSR_TO_COO(CSR A, COO B)

  (* Valeurs identiques *)
  B.nnz := A.nnz
  B.size := A.size
  B.columns := A.columns
  B.values := A.values
  
  (* Remplissage du tableau "rows" *)
  Pour i de 0 à A.size - 1
     Pour k de A.pointer[i] à A.pointer[i+1] - 1
        B.rows[k] := i
     FinPour
  FinPour

FinFonction

complexité: 3*nnz + size = O(nnz + size)

Fonction COO_TO_CSR(COO A, CSR B)
  (* Valeurs identiques *)
  B.nnz := A.nnz
  B.size := A.size
  B.columns := A.columns
  B.values := A.values
  
  (* Remplissage du tableau "pointer" *)
  Si A.nnz > 0
    (* Détermination de B.pointer[i] pour les lignes non vides *)
    i := A.rows[0]
    estVide[i] := FAUX
    B.pointer[i] := 0

    Pour k de 0 à A.nnz - 1
      Si A.rows[k] != i
        i := A.rows[k]
        estVide[i] := FAUX
        B.pointer[i] := k
      FinSi
    FinPour

    (* Remplissage de B.pointer[i] pour les lignes vides *)
    B.pointer[A.size] := A.nnz
    Pour i de A.size - 1 à 0
      Si estVide[i]
        B.pointer[i] := B.pointer[i+1]
      FinSi
    FinPour
  FinSi

FinFonction

complexité: size + nnz + size + size + size + nnz = 2*nnz + 4*size = O(nnz + size)

Fonction COO_TO_CSC(COO A, CSC B)
  (* Valeurs identiques *)
  B.size := A.size
  B.nnz := A.nnz

  (* Determination du nombre d'élements par colonne *)
  Pour j de 0 à A.size - 1
    nb[j] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    nb[A.columns[k]] := nb[A.columns[k]] + 1
  FinPour

  (* Remplissage du tableau "pointer" pour les colonnes non vides *)
  k := 0
  Pour j de 0 à A.size - 1
    Si nb[j] > 0
      B.pointer[j] := k
      k := k + nb[j]
    FinPour
  FinPour

  (* Remplissage du tableau "pointer" pour les colonnes vides *)
  B.pointer[A.size] := A.nnz
  Pour j de A.size - 1 à 0
    Si estVide[j]
      B.pointer[j] := B.pointer[j+1]
    FinSi
  FinPour

  (* Remplissage des tableaux "rows" et "values" *)
  (* nb[j] = nombre d'élément déjà stocké dans la colonne j *)
  Pour j de 0 à A.size - 1
    nb[j] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    j := A.columns[k]
    kB := B.pointer[j] + nb[j]
    B.rows[kB] := A.pointer[k]
    nb[j] := nb[j] + 1
  FinPour

FinFonction

complexité: size + nnz + size + size + nnz = 2*nnz + 3*size = O(nnz + size)
Fonction CSC_TO_COO(CSC A, COO B)
  (* Valeurs identiques *)
  B.size := A.size
  B.nnz := A.nnz

  (* Determination du nombre d'élements par ligne *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    nb[A.rows[k]] := nb[A.rows[k]] + 1
  FinPour

  (* Determination du pointeur sur le premier élément de la ligne i *)
  Si A.size > 0
    start[0] = 0

    Pour i de 0 à A.size - 2
      start[i+1] := start[i] + nb[i]
    FinPour
  FinSi
  
  (* Remplissage des tableaux "columns", "rows" et "values" *)
  (* nb[i] = nombre d'élément déjà stocké dans la ligne i *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour j de 0 à A.size - 1
     Pour k de A.pointer[j] à A.pointer[j+1] - 1
        i := A.rows[k]
        kB := start[i] + nb[i]
        B.rows[kB] := i
        B.columns[kB] := j
        B.values[kB] := A.values[k]
        nb[i] := nb[i] + 1
     FinPour
  FinPour
    
FinFonction

complexité: size + nnz + size + size + nnz = 2*nnz + 3*size = O(nnz + size)
Fonction CSC_TO_CSR(CSC A, CSR B)
  (* Valeurs identiques *)
  B.size := A.size
  B.nnz := A.nnz

  (* Determination du nombre d'élements par ligne *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    nb[A.rows[k]] := nb[A.rows[k]] + 1
  FinPour

  (* Remplissage de "pointer" *)
  B.pointer[0] = 0
  Pour i de 0 à A.size - 1
      A.pointer[i+1] := A.pointer[i] + nb[i]
  FinPour

  (* Remplissage des tableaux "columns" et "values" *)
  (* nb[i] = nombre d'élément déjà stocké dans la ligne i *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour j de 0 à A.size - 1
     Pour k de A.pointer[j] à A.pointer[j+1] - 1
        i := A.rows[k]
        kB := pointer[i] + nb[i]
        B.columns[kB] := j
        B.values[kB] := A.values[k]
        nb[i] := nb[i] + 1
     FinPour
  FinPour
    
FinFonction

Idem pour les fonctions CSC_TO_CSR, CSRtranspose et CSRtranspose.

complexité: size + nnz + size + size + nnz = 2*nnz + 3*size = O(nnz + size)
Fonction transposeCOO(COO A, COO B)
  (* Valeurs identiques *)
  B.size := A.size
  B.nnz := A.nnz

  (* Détermination du nombre d'élements par colonne pour A,
     donc du nombre d'éléments par ligne pour B *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    nb[A.rows[k]] := nb[A.rows[k]] + 1
  FinPour

  (* Détermination du début de la ligne i pour B *)
  Si A.size > 0
    start[0] = 0

    Pour i de 0 à A.size - 2
      start[i+1] := start[i] + nb[i]
    FinPour
  FinSi

  (* Remplissage des tableaux "rows", "columns" et "values" *)
  (* nb[i] = nombre d'élément déjà stocké dans la ligne i de B *)
  Pour i de 0 à A.size - 1
    nb[i] := 0
  FinPour

  Pour k de 0 à A.nnz - 1
    i := A.columns[k]
    kB := start[i] + nb[i]
    B.rows[kB] := i
    B.columns[kB] := A.rows[k]
    B.values[kB] := A.values[k]
    nb[i] := nb[i] + 1
  FinPour
  
FinFonction

B. Algorithmes de répartition

Fonction Scatter(A, nbproc)
  
  nnz = A.nnz
  n = A.size
  nnz_moyen = arrondi(nnz / nbproc)
  rang = 0
  i_fin = 0
  i_debut = 0
  
  TantQue rang < nbproc
    
    Si "Dernier processeur"

      // Envoie toutes les lignes/colonnes restantes
      i_fin := n

    Sinon

      Si
         "Il reste autant de lignes/colonnes que de processeur"
         Ou "La ligne/colonne i_debut contient plus que nnz_moyen éléments"

        // On envoie seulement la ligne i_debut
        i_fin := i_debut + 1

      Sinon

         TantQue
           "Le nombre d'éléments entre les lignes/colonnes i_debut et i_fin - 1
	   est inférieur ou égal à nnz_moyen"
           Et "Il reste au moins une ligne/colonne pour chacun des processeur
	      restants"
           i_fin := i_fin + 1
         FinTantQue

      FinSi

      Envoyer les lignes/colonnes de i_debut à i_fin - 1 au processeur

      i_debut := i_fin
      rang := rang + 1

    FinSi
    

  FinTantQue
  
FinFonction

C. Mesure des performances

figure 12

D. Diagrammes de classes

figure 13

figure 14

Cette page est conforme aux normes du W3C - Auteur : Frédéric WANG - Dernière mise à jour : mercredi 16 septembre 2009