<figure><img src="image_1.png" /></figure>

HAL
open science

# Etude de la parallélisation du produit Matrice/Vecteur creux sur processeurs hétérogènes.

Corentin Rossignon

<figure><img src="image_2.png" /></figure>

## To cite this version:

Corentin Rossignon. Etude de la parallélisation du produit Matrice/Vecteur creux sur processeurs hétérogènes. Calcul parallèle, distribué et partagé [cs.DC]. 2011. hal-00793702

HAL Id: hal-00793702

https://inria.hal.science/hal-00793702v1

Submitted on 31 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Étude sur la parallélisation du produit Matrice/Vecteur creux sur processeurs hétérogènes

Corentin Rossignon

Encadrants INRIA : Olivier Aumage, Samuel Thibault

Encadrant TOTAL : Pascal Hénon

Février - Juin 2011
<table><tr><td>2.9</td><td>Format non testés</td><td>23</td></tr><tr><td>2.9.1</td><td>COO : Coordinate format</td><td>23</td></tr><tr><td>2.9.2</td><td>CSC : Compress Sparse Column</td><td>23</td></tr><tr><td>2.10</td><td>Résumé des performances</td><td>24</td></tr><tr><td><b>3</b></td><td><b>Expérimentation du BCSR</b></td><td><b>25</b></td></tr><tr><td>3.1</td><td>Alignement</td><td>25</td></tr><tr><td>3.2</td><td>Remplissage</td><td>26</td></tr><tr><td>3.3</td><td>Limite du BCSR</td><td>27</td></tr><tr><td>3.4</td><td>Utilisation de GeMV</td><td>28</td></tr><tr><td>3.5</td><td>Changement de granularité</td><td>28</td></tr><tr><td><b>4</b></td><td><b>Conclusion</b></td><td><b>30</b></td></tr></table>
# Introduction

La simulation physique est un élément essentiel de nos jours pour comprendre le monde qui nous entoure. Elle est utilisée aussi bien pour prédire le temps que pour réaliser des essaies nucléaires sans impact radioactif sur l'environnement. Ces simulations sont gourmandes en puissance de calcul et doivent souvent s'exécuter plus vite que dans la nature ; par exemple : mettre deux jours à calculer le temps du lendemain n'est pas intéressant. Les machines réalisant ces calculs doivent donc être de plus en plus puissantes pour pouvoir donner un résultat plus vite et plus précis.

La première solution utilisée jusqu'aux environs de 2005 était d'augmenter la fréquence des processeurs. Cependant, une limite physique a été atteinte autour de 4 Ghz, au-delà la dissipation thermique n'est plus possible. La seconde solution fut d'augmenter le nombre de processeurs. L'inconvénient de cette solution est que le programmeur doit changer sa façon de programmer et souvent il n'arrive pas, ou ne peut pas, utiliser toute la puissance de la machine.

Dernièrement, la programmation GPGPU (General-Purpose Computation on Graphics Hardware) a été proposée, elle consiste à utiliser la carte graphique (Graphics Processing Unit ou GPU) pour effectuer les calculs. Dans certains cas le temps de calcul peut être divisé par 100, mais ici aussi, le programmeur doit programmer différemment avec de nouvelles contraintes mémoires.

Au vu des performances que peut apporter l'utilisation d'un GPU, la société Total, dont les activités se concentrent autour de l'extraction du pétrole brut et de son raffinage, propose d'accélérer leurs simulations en utilisant une solution GPGPU.

Parmi ces simulations, il existe la simulation de réservoir qui est un élément important dans l'étude des champs pétrolifères, elle peut couvrir une période de trente ans. Actuellement cette simulation est faite sur CPU. Ce mémoire propose donc d'étudier l'accélération que l'utilisation hybride GPU - CPU peut apporter par rapport à une solution CPU multi-coeur.
# Chapitre 1

## Contexte

### 1.1 Simulation de réservoir

La simulation de réservoir permet à Total d'optimiser le placement des puits pour l'extraction du pétrole. Dans cette simulation, le temps est discrétisé avec un pas de l'ordre de la journée. Pour chaque pas, il est nécessaire de faire converger un système non-linéaire, pour cela de nombreux systèmes linaires creux doivent être résolus. Les méthodes de résolution de systèmes linaires creux utilisent des matrices.

Les matrices utilisées pour cette simulation sont particulières, elles sont creuses et sont composées de blocs carrés non nuls positionnés sur des diagonales.

<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Matrices utilisées pour la simulation de réservoir.</figcaption></figure>

### 1.2 Les matrices creuses

Les matrices creuses sont des matrices essentiellement composées de zéros, cette propriété peut être utilisée pour réduire son empreinte mémoire ainsi que le nombre d'opérations effectuées lors de l'utilisation d'une telle matrice. Pour le cas du produit matrice/vecteur creux, aussi appelé SpMV, le nombre de multiplication et d'addition est réduit.

Dans cet exemple (voir figure 1.2), seules 4 multiplications et 4 additions sont nécessaires au lieu des 16 multiplications et 16 additions que nous aurions effectuées dans un cas de même taille dense.
<figure><img src="image_5.png" /><figcaption>FIGURE 1.2 – Multiplication matrice/vecteur creux - SpMV.</figcaption></figure>

Pour ne pas avoir à stocker tous les zéros, seuls les éléments non nuls (NZ) sont stockés. Différents formats de stockage de matrices creuses existent comme le COO ou le CSR. Le format a un impact sur la taille de la matrice en mémoire et sur l'ordre des accès mémoires.

## 1.3 Les méthodes itératives

Il existe plusieurs méthodes de résolution de systèmes linéaires creux, dont les méthodes itératives font partie. L'idée principale de ces méthodes est de faire converger la solution en effectuant de nombreux SpMV[Saa96]. Il est aussi possible de préconditionner la matrice pour diminuer le nombre de SpMV à faire.

La méthode ILUk, par exemple, consiste à factoriser partiellement la matrice (A) en deux matrices triangulaires : une inférieure (L) et une supérieure (U). Puis pour résoudre

<figure><img src="image_6.png" /><figcaption>FIGURE 1.3 – Résolution de $Ax = b$ par une méthode ILUk itérative</figcaption></figure>

$Ax = b$, il faut résoudre $Ly = b$ et $Ux = y$ mais comme la matrice n'était que partiellement factorisée, nous devons réitérer les résolutions L et U pour faire converger $x$ vers la solution finale.
## 1.4 La programmation GPGPU

La programmation GPGPU consiste à détourner l'utilisation du processeur graphique (GPU), qui habituellement effectue des calculs liés à l'affichage, pour lui faire effectuer des calculs généralistes. Pour cela, nous devons écrire un code, aussi appelé *Kernel*, dans un langage exprimant le parallélisme (CUDA [NVI11b] pour une carte NVIDIA ou OpenCL [Khr10] pour la portabilité). Dans ces langages, le parallélisme est exprimé sous la forme de Thread, un Thread représente une instance du *Kernel* auquel on a donné un identifiant. Les Threads sont organisés en groupes appelés Warp, les Threads d'un même Warp sont totalement synchrones. Pour CUDA la taille d'un warp est de 32. Les Threads sont répartis dans des blocs (*Block*) et ces blocs sont eux-mêmes répartis dans une grille (*Grid*). La taille des blocs et de la grille déterminent le nombre de Threads. Par exemple, si nous avons 512 Threads par bloc et 512 blocs dans la grille, alors nous aurons un total de 262144 Threads.

<figure><img src="image_7.png" /><figcaption>FIGURE 1.4 – Organisation des Threads pour la programmation GPGPU.</figcaption></figure>

Dans CUDA, l'appel des fonctions GPU et la gestion des identifiants se font de la manière suivante :

```c
__global__ void
exemple(void){
    int identifiant_local = thread Idx.x;
    int identifiant_global = block Idx.x * blockDim.x + thread Idx.x;
}

int
main(){
    exemple<<<512,512 >>>();
    return 0;
}
```

<caption>FIGURE 1.5 – Exemple de code CUDA</caption>

L'identifiant local correspond à l'identifiant d'un Thread à l'intérieur d'un bloc, il est
donc compris entre 0 et 511. Plusieurs Threads dans des blocs différents peuvent avoir le même identifiant local. L'identifiant global est unique à chaque Thread.

<figure><img src="image_8.png" /><figcaption>FIGURE 1.6 – Comparaison architecture CPU - GPU.</figcaption></figure>

Le GPU, à l'inverse du CPU, est surtout composé d'unités de calcul. Cela permet de pouvoir faire s'exécuter un grand nombre de Threads en parallèle. Regardons plus en détail l'architecture des unités de calcul du GPU :

<figure><img src="image_9.png" /><figcaption>FIGURE 1.7 – Architecture d'un Streaming Processor.</figcaption></figure>

Le GPU est composé d'unités SIMD appelées Streaming Processor (SP), ces unités sont capable à un instant T d'exécuter N instructions identiques sur N données potentiellement différentes, avec N la taille d'un warp. En revanche, si un Thread du warp suit un branchement différent des autres, alors il faudra 2 cycles pour exécuter les N instructions.

Il existe deux types de mémoire sur GPU :

- la mémoire partagée (*shared memory*), de petite taille et avec une faible latence partagée entre les Threads d'un bloc.

- la mémoire globale (*global memory*), commune à tous les Threads, mais avec une latence élevée qui peut être masquée par l'utilisation d'un grand nombre de Threads.

Les accès à la mémoire globale présentent une contrainte, pour que les 32 instructions puissent s'exécuter en parallèle, les 32 données doivent être présentes dans le SP, or la mémoire globale ne peut envoyer que 128 octets à la fois (soit 32 float), les Threads doivent donc accéder aux données contiguës en mémoire, ce que l'on nomme faire des
accès coalescés. Si les données sont réutilisées plusieurs fois, il peut être avantageux de les
placer dans la mémoire locale (*shared memory*) pour minimiser la latence et aussi pour
pouvoir faire des accès désordonnés à la mémoire locale plutôt qu'à la mémoire globale.

1.5 StarPU

StarPU[ATNW11] est un support exécutif développé par l'INRIA et le LaBRI. Il fournit un modèle d'exécution unifié pour les machines multi-coeurs équipées d'accélérateurs afin d'exploiter l'intégralité de la puissance de calcul tout en s'affranchissant des difficultés liées à la gestion des données.

L'application doit enregistrer ses données auprès de StarPU qui en interne utilise
une DSM (Distributed Shared Memory) réalisant de manière transparente les transferts
mémoires entre les unités de calcul.

Puis, l'application crée des tâches décrivant pour chaque unité de calcul un code spécifique ainsi que les données nécessaires.

StarPU s'occupera ensuite d'ordonnancer les tâches sur les différentes unités de calcul
en prenant en compte le temps des transferts mémoires et le temps de calcul.

1.6 Solutions existantes

Lors de mes recherches sur les solutions existantes de SpMV sur GPU, deux solutions se sont détachées du groupe : Cusparse et Cusp. Elles utilisent toutes les deux CUDA.

1.6.1 Cusparse

Cusparse[NVI11a] est une bibliothèque développée par NVIDIA et distribuée avec le SDK de CUDA. Cette bibliothèque fournit quelques opérations de bases d'algèbre linéaire creux, tel que le SpMV et le SpMM (produit matrice/matrice creux) pour la catégorie BLAS 2 et BLAS 3. Les BLAS sont des routines d'algèbres linéaires.

En ce qui concerne le SpMV, les matrices sont stockées au format CSR. Plus loin, dans les tests, nous verrons que le format CSR n'est pas le format le plus approprié au SpMV sur GPU.

1.6.2 Cusp

Cusp[BG10][BG08] est une bibliothèque sous licence Apache 2.0, implémentant diverses opérations d'algèbre linéaire creux, dont le SpMV, sur GPU et CPU. Cette bibliothèque est générique (le type des coefficients des matrices sont paramétrables) et performante. Elle propose différents formats de stockage des matrices, tel que COO, CSR, DIA, ELL et un format hybride ELL+COO. Son seul inconvénient est que dans le cas d'une machine multi-GPU, elle ne permet pas d'utiliser tous les GPU en même temps[cus].

Pour notre problème, StarPU doit pouvoir gérer lui-même la mémoire. Or, Cusp ne nous laisse pas la gestion des transferts mémoires en plus de ne pas être multi-GPU ; il n’est donc pas possible d’utiliser cette bibliothèque en même temps que StarPU. Malgré cela
Cusp est performante dans le cas de matrices génériques. C'est pourquoi nous l'utiliserons comme référence de performances pour nos *Kernels*.

## 1.7 Proposition

Notre objectif est de pouvoir faire un produit matrice/vecteur creux, optimisé pour nos matrices, qui utilise à la fois les unités de calculs CPU et GPU. Pour cela nous allons découper nos matrices en plusieurs morceaux représentant chacun une tâche StarPU, ces tâches pourront être ensuite réparties sur les différents cœurs de calcul.

Pour les unités GPU, nous rechercherons un format adapté à nos matrices qui pourra être moins générique mais plus performant que les formats utilisés par Cusp.
# Chapitre 2

## Formats de stockage et
performances

Le format de stockage des matrices creuses ayant un impact sur les performances du
SpMV, nous avons essayé de déterminer quel format est le plus performant sur GPU pour
nos matrices.

### 2.1 Réalisation des tests

La machine de test est composée de 3 Tesla C2050 (puissance de crête : 1030 Gflops et 3 Go de mémoire pour un GPU) et de 2 Xeon X5650 (6 cœurs, 2.66 GHz pour un CPU).

Les résultats sont exprimés en flops, pour chaque élément non nuls de la matrice, nous effectuons une multiplication et une addition, la formule pour calculer les flops est donc :
$$ (2 * NNZ) / Temps, $$
NNZ représente le nombre d'éléments non nuls.

Les matrices utilisées sont générées par un programme développé par Total et sont représentatives de ce qui est réellement utilisé.

Ces matrices représentent un domaine 3D. Pour chaque élément du domaine, le nombre de variables de l'élément correspond à la taille des blocs de la matrice. La dimension de la matrice est égale à la dimension du domaine à la puissance 3.

Dans les résultats, nous nommerons les matrices de la façon suivantes : "NNZ TailleBloc TailleDimension".

La taille des matrices choisies correspond à ce qui est actuellement utilisé dans la simulation de réservoir. Dans les cas réels, le paramètre TailleBloc varie entre 3 et 21, nous avons choisi les valeurs 4, 8, 9 et 16 pour pouvoir évaluer l'impact de la taille des blocs sur les performances de chaque format. La valeur 9 est particulière car il s'agit d'une valeur qui dans le cas du BCSR sur GPU peut donner de mauvaises performances. La taille des dimensions ne fait qu'augmenter la taille de la matrice, en prenant 20 et 30 nous aurons à la fois des problèmes de petite taille et des problèmes de taille raisonnable. Nous avons aussi rajouté une matrice composée de blocs carrés 16x16 et 40 en dimension, ainsi nous aurons une matrice de grande taille qui nous permettra d'estimer les limites des formats.
## 2.2 SpMV sur CPU

Pour pouvoir avoir une référence, nous utilisons la bibliothèque MKL d'Intel regroupant des routines d'algèbres linéaires efficaces. Nous utiliserons la fonction MKL qui effectue un produit matrice/vecteur creux sur une matrice stockée au format CSR.

<figure><img src="image_10.png" /><figcaption>FIGURE 2.1 – Résultats du SpMV sur 2 Xeon X5650 de 6 cœurs chacun soit 12 cœurs cadencés à 2.66 GHz.</figcaption></figure>

Si nous utilisons un coeur de calcul CPU nous obtenons environ 1 Gflop. Dans le meilleur cas, avec 12 cœurs de calcul nous pourrions espérer avoir 12 Gflops. Cependant, le débit mémoire nous limite à 25 Go/s (débit théorique : 32 Go/s), ce qui correspond à 5.5 Gflops. Par conséquent, augmenter le nombre de CPU serait inutile.

Nous venons de voir qu'au maximum nous obtenons 5.5 Gflops, ce qui reste peu par rapport au 15 Gflops que nous pourrions obtenir sur un seul GPU[VOFG10]. Nous décidons donc que pour notre problème, le CPU n'est pas assez performant pour aider le GPU. A présent, nous travaillerons uniquement sur GPU.

## 2.3 SpMV sur GPU

Pour l'implémentation d'un *Kernel* SpMV sur GPU, plusieurs questions doivent être posées :

- Combien de Threads doivent être créer ?

- Comment distribuer le travail entre les Threads ?

- Les accès mémoires sont-ils coalescés ?
Pour chaque format nous essayerons de trouver la meilleure réponse.

Pour réaliser le SpMV sur plusieurs unités de calcul, nous essayons différentes granularités de découpage des matrices :

- 1 partie : dans ce cas la matrice n'est pas découpée, cela nous permet de mesurer l'efficacité du SpMV sur un seul GPU.

- 3 parties : nous pouvons maintenant distribuer les calculs sur chacun de 3 GPUs et mesurer si le format peut être découpé efficacement.

- 9 parties : nous créons un peu plus de tâche que de GPU, cela nous permet de voir si le surcoût du découpage n'est pas trop important pour pouvoir avoir des tâches plus petites pour que StarPU puisse les ordonnancer efficacement sur des GPUs hétérogènes.

- 30 parties : nous avons 10 fois plus de tâches que de GPU, nous pouvons voir si le format fonctionne aussi bien sur un petit ensemble de données pour les mêmes raisons que 9 parties.

Les tests suivants ont été réalisés en simple précision.

## 2.4 Compressed Sparse Row

### 2.4.1 Analyse du format

<figure><img src="image_11.png" /><figcaption>FIGURE 2.2 – Compressed Sparse Row.</figcaption></figure>

Les éléments non nuls sont stockés par ligne (NZ), pour chaque élément nous avons sa colonne (COL) et sa ligne doit être calculée en utilisant PTR. Le tableau PTR correspond à l'indice dans NZ du premier élément d'une ligne.

### 2.4.2 Implémentation sur GPU

Les Threads peuvent être distribués de deux manières : 1 Thread par ligne ou 1 warp (32 Threads) par ligne. La solution d'un Thread par ligne ne permet pas d'avoir des accès mémoires coalescés à NZ et COL, et donc produit de très mauvaises performances. A l'inverse, la solution du warp par ligne permet d'avoir des accès mémoires coalescés à NZ et COL, mais oblige une réduction du résultat des 32 Threads à la fin du calcul qui peut être fait efficacement avec la technique de Mark Harris[Har]. Comme l'algorithme fonctionne
par ligne, la matrice peut être découpée en nombre égal de lignes. Le vecteur multiplicatif doit être transmis à chaque unité de calcul et le vecteur résultat est découpé comme la matrice. Quelque soit la méthode utilisée, l'ordre des accès au vecteur multiplicatif reste indéfini.

<figure><img src="image_12.png" /><figcaption>FIGURE 2.3 – Répartition des Threads dans le format CSR.</figcaption></figure>

### 2.4.3 Résultats

<figure><img src="image_13.png" /><figcaption>FIGURE 2.4 – Résultats du SpMV pour le format CSR.</figcaption></figure>
Nous pouvons remarquer que la taille des blocs de la matrice est un facteur important, cela vient des accès mémoires, une taille de blocs se rapprochant d'un multiple de 32 permet d'augmenter le nombre d'accès coalescés au vecteur multiplicatif. Les blocs ayant un impact sur le format, nous essayerons plus tard le format BCSR qui utilise des blocs plutôt que des coefficients.

Le nombre de découpages n'a pas un fort impact sur les performances car il n'y a aucune dépendance entre les morceaux du vecteur résultat. Il suffit juste de recoller les bouts. Les résultats mono GPU sont aussi bons qu'avec Cusp. Nous obtenons un speed up de 2.8 en utilisant 3 GPUs, ce qui montre que ce format se découpe bien.

## 2.5 Diagonal format

### 2.5.1 Analyse du format

<figure><img src="image_14.png" /><figcaption>FIGURE 2.5 – Diagonal format.</figcaption></figure>

Le format DIA tient compte de la structure de la matrice, si les éléments non nuls de la matrice sont concentrés sur des diagonales, ce qui est le cas dans les matrices que nous voulons manipuler, alors il suffit de stocker toutes les diagonales qui contiennent un élément non nul en remplissant le vide par des zéros et dans un autre tableau, stocker le numéro de la diagonale, les nombres négatifs correspondent aux diagonales basses, les positifs aux diagonales hautes, le 0 à la diagonale centrale. Grâce au numéro de diagonale, nous pouvons retrouver la taille de la diagonale ($nx - |diag|$).

### 2.5.2 Implémentation sur GPU

Comme pour le CSR nous avons plusieurs choix de distribution ; 4 au total, 1 Thread/-warp par ligne/diagonale. Toutes les solutions de distribution par diagonale impliquent soit d'avoir un vecteur résultat temporaire par Thread/warp soit d'effectuer des écritures atomiques, ces deux solutions sont peu performantes voir impossibles. En attribuant un Thread à une ligne, et grâce à la synchronisation du warp, nous sommes sûrs que tous les Threads d'un warp travaillent sur la même diagonale et accèdent donc de manière coalescer et directe à NZ et au vecteur multiplicatif.
<figure><img src="image_15.png" /><figcaption>FIGURE 2.6 – Répartition des Threads dans le format DIA.</figcaption></figure>

Toutefois, cette solution comporte un problème sur les bords de la matrice, il y aura des Threads qui n'auront rien à faire, par exemple pour la diagonale -31, seul le Thread 32 aura une opération à effectuer, les 31 autres Threads (1-31) attendront.

### 2.5.3 Résultats

<figure><img src="image_16.png" /><figcaption>FIGURE 2.7 – Résultats du SpMV pour le format DIA classique.</figcaption></figure>

Si nous comparons nos résultats avec ceux de Cusp, qui sont environ 10 Gflops supérieurs, nous pouvons voir que nous n'obtenons pas un résultat optimal. Le problème provient du stockage de la matrice. En effet, au lieu de stocker des diagonales avec différentes tailles,
nous devrions normaliser leurs tailles pour permettre au GPU de trouver plus facilement
une diagonale (évite l'utilisation d'un tableau d'indirection). De plus, la normalisation
des diagonales n'augmente que la taille en mémoire de la matrice sans changer le nombre
d'opérations à effectuer.

<figure><img src="image_17.png" /><figcaption>FIGURE 2.8 – Résultats du SpMV pour le format DIA orienté GPU.</figcaption></figure>

Les résultats sont meilleurs et concordent avec Cusp. Cependant la structure de nos
matrices n’est pas optimale pour cette version du DIA, cela provient du fait qu’elles sont
blocs-diagonales. Le problème apparaît pour les diagonales du bord des blocs, car ces
diagonales se retrouvent avec plus de zéro que d’éléments non nuls.

<figure><img src="image_18.png" /><figcaption>FIGURE 2.9 – Zoom sur le problème de remplissage d'une diagonale.</figcaption></figure>

Pour résoudre ce problème, il faudrait créer un nouveau format qui serait un DIA par
bloc.
## 2.6 ELLPACK format

### 2.6.1 Analyse du format

<figure><img src="image_19.png" /><figcaption>FIGURE 2.10 – ELLPACK format.</figcaption></figure>

Le ELL consiste d'abord à "pousser" tous les éléments vers la gauche, puis à remplir de zéros pour créer une structure rectangulaire. Pour chaque élément nous enregistrons sa colonne dans COL (ou -1 pour les zéros). Ce format peut être stocké soit au format *row major* (les données contiguës appartiennent à la même ligne), soit au format *column major* (les données contiguës appartiennent à la même colonne).

### 2.6.2 Implémentation sur GPU

La première question à se poser concerne l'ordre de stockage de la matrice. Le *row major* revient à avoir un CSR avec des zéros supplémentaires, alors que le *column major* nous offre plus de parallélisme en permettant d'avoir un Thread par ligne tout en gardant des accès coalescés à NZ.

<figure><img src="image_20.png" /><figcaption>FIGURE 2.11 – Répartition des Threads dans le format ELL.</figcaption></figure>
### 2.6.3 Résultats

<figure><img src="image_21.png" /><figcaption>FIGURE 2.12 – Résultats du SpMV pour le format ELL en column major.</figcaption></figure>

Les résultats sont plutôt bons car nous travaillons sur des matrices régulières, il n'y a donc quasiment pas d'ajout de zéros. De plus, grâce à la structure en bloc de nos matrices, les accès au vecteur multiplicatif se retrouvent mutualisés. Les performances sont aussi bonnes que Cusp.

## 2.7 Jagged Diagonal Storage

### 2.7.1 Analyse du format

Le JAD est un format dérivé du ELL, pour éviter l'ajout de zéros les lignes sont réordonnées par taille et les permutations sont stockées dans PERM. Le vecteur final s'obtient en appliquant PERM au vecteur résultat.

### 2.7.2 Implémentation sur GPU

Le JAD étant un dérivé du ELL, les problématiques de l'implémentation sont les mêmes que pour le ELL, nous avons donc un Thread par ligne.

### 2.7.3 Résultats

Les résultats (voir figure 2.15) sont aussi bons que ceux du format ELL. Cependant, le JAD peut se montrer plus performant que le ELL dans le cas de matrices non régulières,
mais ce n’est pas notre cas.

Cusp n’implémente pas le format JAD. En contrepartie, un format hybride ELL+COO est utilisé pour les matrices non régulières.

<figure><img src="image_22.png" /><figcaption>FIGURE 2.13 – JAgged Diagonale storage.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 2.14 – Répartition des Threads dans le format JAD.</figcaption></figure>

## 2.8 Block Compress Sparse Row

### 2.8.1 Analyse du format

Le BCSR reprend le principe du CSR mais stocke des blocs. Pour minimiser l’ajout de zéros, nous avons des blocs à une dimension, ce qui correspond à regrouper des lignes. A l’intérieur des blocs, les données peuvent être stockées soit en *column major* soit en *row major*.
<figure><img src="image_24.png" /><figcaption>FIGURE 2.15 – Résultats du SpMV pour le format JAD.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 2.16 – Block Compress Sparse Row. Taille des blocs = 2.</figcaption></figure>

## 2.8.2 Implémentation sur GPU

Le premier choix d'implémentation consiste à choisir un Thread par bloc ou un warp par bloc. Le choix du warp permet d'exploiter le parallélisme offert par les blocs, mais impose une limitation. Pour que les Threads n'effectuent qu'une seule réduction, ils doivent rester sur la même ligne durant le calcul d'un bloc, ce qui oblige à avoir un diviseur de 32 comme taille de bloc.

Dans le cas où le BCSR est constitué de blocs dont la taille n'est pas un diviseur de 32, par exemple 3, la répartition du warp dans le bloc pose problème et oblige l'utilisation
<figure><img src="image_26.png" /><figcaption>FIGURE 2.17 – Répartition des Threads dans le format BCSR.</figcaption></figure>

pour chaque thread d'un tableau temporaire sur lequel une réduction sera faite à la fin.

Ensuite nous devons choisir la taille des blocs, celle-ci dépendant de la taille des blocs en entrée et ne pouvant pas se calculer, nous essayerons dans le Chapitre 3 différentes tailles de blocs.

### 2.8.3 Résultats

Les résultats sur les figures 2.18 et 2.19 ont été obtenus avec des blocs BCSR de taille 8.

L'organisation des données *column major* offre de meilleures performances qu'en ligne, cela s'explique par le fait que les accès au vecteur multiplicatif deviennent mutualisés tout en gardant des accès à NZ et COL coalescés. Nous pouvons remarquer que la taille des blocs de la matrice a un impact important sur les performances. Dans ce cas, il faut avoir une taille de blocs multiple de 8 pour obtenir le maximum de performances. Malgré les mauvaises performances obtenues avec les blocs 9x9, le BCSR reste très efficace.
<figure><img src="image_27.png" /><figcaption>FIGURE 2.18 – Résultats du SpMV pour le format BCSR en *row major*.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 2.19 – Résultats du SpMV pour le format BCSR en *column major*.</figcaption></figure>
## 2.9 Format non testés

Les formats suivants n'ont pas été testé car leurs implémentations sur GPU seraient inefficaces.

### 2.9.1 Coordinate format

<figure><img src="image_29.png" /><figcaption>FIGURE 2.20 - Coordinate format.</figcaption></figure>

Ce format stocke, pour chaque élément non nul (NZ) de la matrice, ses coordonnées (ROW, COL) ainsi que sa valeur. L'ordre des éléments n'étant pas spécifié, il est impossible de savoir où les Threads vont écrire. Il faudrait donc créer un vecteur résultat par Thread puis faire une réduction sur tous ces vecteurs (utilisation coûteuse de la mémoire) ou bien effectuer des opérations atomiques sur le vecteur résultat (sérialisation du code).

### 2.9.2 Compress Sparse Column

<figure><img src="image_30.png" /><figcaption>FIGURE 2.21 - Compress Sparse Column.</figcaption></figure>

Pour le CSC, les valeurs sont stockées par colonne, à l'inverse du CSR. Le problème du portage sur GPU est très proche du problème rencontré pour le COO. Ici, il faudrait 1 vecteur par warp (encore trop coûteux) ou alors utiliser des opérations atomiques.
## 2.10 Résumé des performances

Le CSR est un format efficace sur CPU mais n'est que relativement efficace sur le GPU.

Le DIA doit être un peu modifié pour le GPU, dans le cas d'une matrice composée uniquement de diagonales simples, il se révélera être le plus efficace, mais dans notre cas, les blocs induisent un surplus de calcul. Un DIA par bloc pourrait résoudre notre problème.

Le ELL et le JAD sont aussi efficaces sur GPU mais le JAD est à privilégier dans le cas d'une matrice non régulière.

Pour finir, le BCSR *column major* se révèle jusqu'à 1,5x plus performant que les autres formats, et même dans le cas d'une matrice avec des blocs 9x9 il reste le plus performant.
# Chapitre 3

# Expérimentation du BCSR

L'algorithme actuel du BCSR nécessite des blocs dont la taille est une puissance de 2 pour pouvoir être efficace. Or, dans la cas des matrices que nous étudions, nous n'aurons pas une puissance de 2 en taille de blocs. Nous proposons d'étudier le comportement de l'algorithme dans ces cas là. Pour cela nous avons généré des matrices dont la taille des blocs est comprise entre 2 et 32. Pour chaque matrice nous avons choisi un paramètre dimension de telle sorte que le nombre d'élément non nuls soit d'environ 27 Millions.

Pour chaque matrice nous avons testé deux versions du BCSR en faisant varier la taille interne des blocs BCSR parmi 2, 4, 8 et 16.

## 3.1 Alignement

Si la taille des blocs de la matrice en entrée n'est pas un multiple de la taille des blocs BCSR choisi, alors l'algorithme fait du remplissage de zéros pour obtenir un bloc BCSR.

<figure><img src="image_31.png" /><figcaption>FIGURE 3.1 – Exemple d'un défaut d'alignement pour une matrice bloc 3 et un BCSR bloc 4</figcaption></figure>

Dans le cas d'une matrice composée de blocs carrés de taille 3 et d'un BCSR avec une taille de bloc 4, nous avons un défaut d'alignement pour chaque bloc ce qui conduit à presque multiplier par 2 la taille de la matrice.

Lorsque les blocs sont alignés, il n'y a aucun remplissage, par contre dans les autres cas il y a remplissage, ce qui conduit à augmenter le nombre de calcul.
<figure><img src="image_32.png" /><figcaption>FIGURE 3.2 – SpMV - BCSR col - 3 GPU - Alignement - 27 Mnnz ± 1 Mnnz.</figcaption></figure>

Pour un BCSR avec des blocs de taille 8 nous avons un cas particulier, si la matrice a des blocs de taille 12 alors le remplissage ne se fera qu'une fois sur deux, réduisant ainsi le nombre d'opération.

Quelque soit le nombre de blocs de la matrice, l'utilisation du BCSR en bloc 2 donne toujours les plus mauvaises performances alors que le BCSR en bloc 8 représente le meilleur compromis.

## 3.2 Remplissage

Pour limiter le remplissage, nous rajoutons des lignes de zéros dans chaque bloc pour que celui-ci obtienne un nombre de lignes en puissance de 2. Cela permet de supprimer le remplissage en $x$ que nous avions.

L'ajout de lignes agrandit le vecteur résultat. A la fin du calcul du SpMV, nous devrons donc retirer les résultats inutiles du vecteur.

Dans cet exemple notre matrice blocs carrés est transformée en matrice avec des blocs 3x4. Nous avons donc minimisé l'ajout de zéros.

Les résultats sont plutôt bons lorsque la taille des blocs de la matrice est légèrement inférieure à la puissance de 2 supérieure. Par contre si la taille est légèrement supérieure il vaut mieux utiliser la méthode du défaut d'alignement.
<figure><img src="image_33.png" /><figcaption>FIGURE 3.3 – Exemple d’un remplissage pour une matrice bloc 3 et un BCSR bloc 4.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 3.4 – SpMV - BCSR col - 3 GPU - Remplissage - 27 Mnnz ± 1 Mnnz.</figcaption></figure>

## 3.3 Limite du BCSR

Nous cherchons maintenant les limites du format BCSR, pour cela nous choisissons une matrice composée de blocs carrés 16x16 et nous faisons varier sa dimension entre 1 et 50.

Nous obtenons une limite autour de 140 Gflops lorsque nous utilisons 3 GPUs. Dans les tests que nous avons effectués, le choix d'une matrice de dimension 40 nous donne presque les meilleures performances.
<figure><img src="image_35.png" /><figcaption>FIGURE 3.5 – SpMV - BCSR col - 3 GPU - Matrice bloc 16 x 16.</figcaption></figure>

## 3.4 Utilisation de GeMV

Sur CPU, lorsqu'une matrice creuse comporte des parties denses assez grandes, il peut être avantageux d'utiliser une multiplication matrice/vecteur dense, ou GeMV, sur les parties denses. Le GeMV est un BLAS 2.

Dans notre cas pour créer des blocs denses, nous devons regrouper les blocs d'une même ligne. Il nous faut ensuite construire le vecteur à multiplier pour chaque bloc dense.

Dans le cas où l'on doit exécuter plusieurs SpMV à la suite sur une même matrice, la compression de la matrice n'est faite qu'une seule fois mais il faut créer à chaque itération autant de vecteurs que de GeMV à effectuer. Dans le cas du GPU, il faut envoyer ces vecteurs sur le GPU, et faire un appel GeMV pour chacun. Or si nous souhaitons utiliser le GeMV fourni par NVIDIA dans Cublas, il faut faire les transferts et les appels à GeMV un à un.

L'efficacité du calcul ne parvient pas à amortir le coût du transfert et des appels.

## 3.5 Changement de granularité

Nous avons aussi essayé de modifier le nombre de Threads affectés à un bloc BCSR. Dans un premier temps, nous avons mis 64 Threads par bloc, pour cela nous avons été obligé de rajouter des barrières de synchronisation aux endroits où l'on supposait auparavant les Threads totalement synchrones. Nous obtenons une performance de 70 Gflops contre 100 Gflops pour le même problème avec 32 Threads.
Dans un second temps, nous avons testé l'inverse, nous avons attribués 16 Threads par bloc, le code n'a pas à être changé, mais les accès mémoire ne sont plus coalescés. Les performances sont de 83 Gflops.
# Chapitre 4

## Conclusion

Le produit matrice/vecteur creux présente de nombreuses difficultés pour être porté sur GPU. De part son irrégularité, nous sommes obligés d'utiliser des tableaux d'indirections pour accéder aux éléments. Or l'utilisation de ces tableaux sur GPU est très coûteuse, c'est pour cela que nous cherchons un format adapté aux GPUs.

Notre recherche d'un format adapté aux matrices blocs-diagonales, nous a conduite au format BCSR dans lequel les blocs sont rangés par colonne. Malgré les bonnes performances du BCSR, les performances de crête du GPU sont loin d'être atteintes. Cela s'explique par le fait que nous faisons trois accès mémoire, dont un indirect, pour seulement deux opérations. Un format DIA par bloc pourrait supprimer cet accès indirect et donc améliorer les performances, il sera prochainement testé.

L'utilisation du CPU pour aider le GPU s'est révélée contre performante, le CPU n'effectue pas assez vite la multiplication pour couvrir le coût des transferts mémoires. Pour pouvoir accélérer la simulation de réservoir, nous devrons porter sur GPU toute la méthode itérative de résolution de systèmes linéaires creux. Cela nous permettra de transférer la matrice qu'une seule fois au début de chaque résolution.

Parmi les formats *standards* que nous avons essayés, le ELL et le JAD se sont montrés efficaces mais moins que le BCSR et en plus oblige une réduction du vecteur résultat là où le BCSR ne fait que rassembler les morceaux.

Dans la poursuite de mon stage, nous nous concentrerons sur l'ajout du format BCSR dans le code de Total, puis nous implémenterons le reste de la résolution de systèmes linéaires creux sur GPU[LS10] pour éviter les transferts mémoires $CPU \Leftrightarrow GPU$.
Bibliographie

[ATNW11] Cédric Augonnet, Samuel Thibault, Raymond Namyst, and Pierre-André Wacrenier. StarPU : A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures. Concurrency and Computation : Practice and Experience, Special Issue : Euro-Par 2009, 23 :187–198, February 2011.

[BG08] Nathan Bell and Michael Garland. Efficient sparse matrix-vector multiplication on cuda. December 2008.

[BG10] Nathan Bell and Michael Garland. Cusp : Generic parallel algorithms for sparse matrix and graph computations, 2010. Version 0.1.2; http://cusp-library.googlecode.com Online : accessed 11 juin 2011.

[cus] Issue 60 : Does cusp automatically support multiple gpu. http://code.google.com/p/cusp-library/issues/detail?id=60 Online : accessed 11 juin 2011.

[Har] NVIDIA Developer Technology Mark Harris. Optimizing parallel reduction in cuda.

[Khr10] Khronos OpenCL Working Group. The OpenCL Specification, version 1.1.36, September 2010. www.khronos.org/registry/cl/specs/opencl-1.1.pdf Online : accessed 11 juin 2011.

[LS10] R. Li and Y. Saad. GPU-accelerated preconditioned iterative linear solvers. Technical Report umsi-2010-112, 2010.

[NVI11a] NVIDIA. CUSPARSE User Guide, version 4.0, January 2011. http://developer.download.nvidia.com/compute/DevZone/docs/html/CUDALibraries/doc/CUSPARSE_Library.pdf Online : accessed 11 juin 2011.

[NVI11b] NVIDIA. NVIDIA CUDA C Programming Guide, version 4.0, May 2011. http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_Programming_Guide.pdf Online : accessed 11 juin 2011.

[Saa96] Yousef Saad. Iterative Methods for Sparse Linear Systems. PWS, Boston, 1996.

[VOFG10] F. Vazquez, G. Ortega, José-Jesús Fernández, and Ester M. Garzón”. Improving the performance of the sparse matrix vector product with GPUs. IEEE Computer Society, 2010.