Les recherches dans l'accélération des calculs de mécaniques des fluides nous ont conduit à considérer le transfert d'une partie de ces codes sur des processeurs graphiques. On peut désormais les voir comme des machines vectorielles massivement parallèles.
Ce domaine a été très actif depuis deux ans et l'évolution technologique a poursuivi son rythme effréné. Nous avons mené un travail de bibliographie, de veille technologique et d'instrumentation pour pour pouvoir dresser un paysage complet.
Nous présentons, dans la première partie, l'architecture des cartes graphiques unifiées. En deuxième partie, les outils de programmation disponibles. Les performances et les précision de la multiplication de matrices sont aux troisième chapitre. En quatrième partie nous résolvons les équations de Laplace, se plaçant dans une optique de scalabilité. De futurs axes de recherches sont évoqués en dernier lieu.
La simulation numérique
La monté en puissance des ordinateurs a introduit une rétroaction nouvelle de la technique informationnelle sur de vastes domaines des techniques de la matière et de l'énergie. La simulation numérique peut être définie comme une démarche qui consiste à construire, dans la mémoire d'un ordinateur, une image numérique d'un système physique qui peut être par exemple, un volume de fluide. Cette image comporte schématiquement deux éléments : d'une part les paramètres physiques, pression, température, vitesse, qui en chaque point définissent l'état du fluide à un instant donné - ce sont les conditions initiales ; d'autre part l'ensemble des équations qui régissent l'évolution du fluide. Il s'agit dans ce cas, d'un ensemble d'équation différentielles non linéaires, les équations de Navier Stokes. L'ordinateur est programmé de façon à faire évoluer l'image de l'état du fluide en respectant les contraintes imposées par les équations qui régissent cette évolution.Quel sont les usages de cette technique ? En premier lieu elle permet de prévoir l'évolution d'un système physique en l'assimilant à l'évolution de son image numérique. Rien ne contraint à faire évoluer l'image au même rythme que l'objet, et une évolution accélérée crée la possibilité de prévoir. C'est ce que l'on fait dans la prévision météorologique où l'objet est l'atmosphère terrestre.
La mécanique des fluides offre un domaine de prédilection pour l'usage de cette méthode : si les équations qui gouvernent le comportement des fluides sont parfaitement connues, elles ne se prêtent pas à accéder directement à la compréhension des mouvements du fluide. Le détour par la simulation numérique est le seul accès viable à cette connaissance [1].
La science moderne est faite d'observation, de théories, d'expériences... et de simulations numériques.
Motivations
Dans les années 2000, la tendance affiché par les prédictions de la loi de Moore s'est infléchie. Si les projections s'était poursuivies au rythme qu'elle connurent dans les années 90, nous serions équipés en 2010 de processeurs cadencés à 10 GHz. Or la réduction de la taille des transistors et l'augmentation des fréquence d'horloge atteint ces limites, dues aux phénomènes de transferts et de dissipations caloriques. Aux échelle nanométriques, des problèmes de nature électronique se présentent, qui poussent à explorer d'autre voies. Le processeur Von Neumann n'a pas été conçu pour suivre l'évolution de l'espace disponible.
L'autre difficulté surgit dans l'excès de consommation électrique des puces : un processeur dont le cycle d'horloge oscille à un gigahertz réclame environ 25 watts pour fonctionner, tandis qu'environ 100 watts sont nécessaires pour faire fonctionner un processeur cadencé à 3 GHz.
L'étendue du parc informatique mondial nous amène à reconsidérer le bilan énergétique produit par les PC, les serveurs webs, les DNS, les moteurs de recherches et les centres de calculs. Selon Wire Magazine, rapportant une étude d'AMD [2], l'électricité totale consommée par les moteurs de recherches les plus important en 2006 avoisine 5 gigawatts.
Puisque les CPUs ne peuvent être accélérés plus avant, les efforts se sont tournés sur le parallélisme des unités de traitements. C'est actuellement le seul moyen dont nous disposons pour mettre à contribution davantage de puissance de calcul [3]. La distribution des calculs sur des ordinateurs massivement parallèle, l'étude des machines et des algorithmes parallèle n'est plus un exercice académique mais une nécessité.
En l'état actuel, l'algèbre linéaire se manifeste comme un candidat idéal à l'accélération de code de calculs à grands nombres d'équations. Si le hardware évolue rapidement, mute, parfois devient obsolète, le software, lui, reste et poursuit la logique pour laquelle il a été conçu. Comme on ne peut pas réécrire tous les codes pièces par pièces, on se tourne vers des solutions d'accélération intermédiaire : en remplaçant le matériel existant et les librairies logicielles par des contreparties fonctionnant sur des architectures hybrides, mais tout en conservant les codes historiques. A terme, il sera possible de mettre à la disposition des chercheurs et des ingénieurs, l'équivalent de supercalculateur dans leur station de travail.
Ce phénomène n'est pas nouveau, l'architecture des supercalculateurs inspire celle des microprocesseurs : les techniques mises en oeuvre dans les microprocesseurs sont en premier lieu des principes macroscopiques. En ce sens, Les théories de la programmation parallèles sont plus vérifiables que jamais. Les techniques de calcul parallèle existantes, dans des outils tels que MPI ou OpenMP, à l'échelle des clusters pour le premier et des machines multiprocesseurs pour le second, ont été adaptées, moyennant quelques contraintes, aux GPU. En bref, nous insistons sur un fait : le matériel évolue, mais le logiciel existe pour toujours. L'intrication logicielle/matérielle, telle qu'elle existe dans les systèmes d'exploitation par exemple, est levée dans les couches logicielles relevant de la logique métier, mais revient systématiquement dès lors qu'il sera question d'optimisation et d'arithmétique des ordinateurs, ce qui ne manque par d'arriver dans les codes de calculs.
Les processeurs graphiques
Les unités de traitement graphiques (Graphical Processing Units, GPU) sont à l'origine dédiés aux affichages et aux rendus de géométries et d'images. Leader dans ce domaine, l'industrie du jeu vidéo à guidé leur évolution. A ce titre, depuis 1992, le comité de l'Architecture Review Board (ARB), réunit un consortium des principaux fabriquants de matériels.
Pour que le développement logiciel puisse reposer sur des standards l'ARB s'est entendu afin d'unifier le modèle de programmation des cartes graphiques, par un versionnage et des test de compatibilités rigoureux de ses standards.
Parmis ses résultats, la spécification d'OpenGL a apporté des solutions d'émulations logicielles, qui, au fil des versions, ont été intégrées aux matériels. Une parité logicielle/matérielle a émergé, l'amélioration du matériel permettant la mise au point de nouvelles solutions logicielles, solutions qui a leur tour s'intègrent au matériel.
Les GPU, présents dans chaque ordinateur, ont évolués au tournant de la décennie en coprocesseurs d'arithmétiques vectorielles massivement parallèles, ce qui explique aujourd'hui leur présence dans les domaines du calcul scientifique, numéro un du parallélisme de masse.
Enjeux industriels
Délivré en 2006 après une longue période de gestation en bêta, CUDA a conquis de nombreux territoires dans le domaine de l'ingénierie scientifique. Depuis, Nvidia a considérablement revu la conception et le positionnement de ses GPU et la désignation commune doit les reconnaître comme des appareils versatiles dont l'usage dépasse celui des jeux électroniques et du graphisme.
Pour Nvidia, le HPC est à la fois une opportunité de vendre plus de puces et une assurance contre le futur incertain des GPU discrets. Bien que les les cartes graphiques de Nvidia aient longtemps été prisées par les joueurs, le marché a changé. Lorsque AMD acquis ATI en 2006, Nvidia demeura le plus gros vendeur indépendant de GPU. En réalité le "géant vert" est même le seul vendeur indépendant, les autres compétiteurs ayant échoué ou été racheté au cours des ans. Le statut de dernier survivant est enviable à condition que le marché pour les GPU discrets demeure stable [4]. AMD et Intel projettent l'intégration de cores graphiques dans leur processeurs futurs ce qui pourrait considérablement réduire le marché des GPU discrets. Annoncé pour la fin 2009, le projet Larabee d'Intel a été repoussé à une date inconnue. Les différents échanges commerciaux qui ont suivis laissent supposer qu'Intel ne sortira pas de cartes graphiques externes.
L'intérêt de Nvidia pour développer sa recherche est manifeste : une tactique de collecte des "success story" [5,6] dans l'accélération des codes grâce à sa technologie fait non seulement une bonne publicité, certes, la place aussi dans une situation favorable pour déterminer les projets profitables au titre du commerce par l'adaptation au marché, mais aussi et surtout par l'amélioration de sa propre R&D.
Un exemple : le GPU permet d'accélérer les codes de simulations numériques de dissipations thermiques, qui a leur tour ouvre une nouvelle perspectives dans l'amélioration des chipsets. Le temps de prototyper et la boucle de rétroaction logicielle/matérielle se referme en quelques semaines pour le constructeur, en plusieurs mois pour ses clients, et l'écart entre les deux demeure suffisement important pour que la logique commerciale s'intègre dans le paysage de la stratégie américaine consistant à maintenir cet écart technologique avec le reste du monde.
Apple est l'un des premiers à proposer une gamme d'ordinateurs équipés de deux processeurs graphiques, l'un pour l'affichage, l'autre pour l'accélération d'OSX et de ses applications. Apple supporte vivement la spécification du standard OpenCL.
Microsoft, sorti de l'ARB en 2003, a fait chambre à part en poursuivant l'évolution de DirectX11. Sans grande surprise, les innovations sont rares et inspirées de travaux existants. Windows demeurant la plateforme de prédilection de l'industrie du jeu vidéo, les pilotes graphiques vont poursuivre encore longtemps leur interface avec le système d'exploitation de Redmond.
On peut supposer que les GPU seront neuf fois plus rapides dans trois ans.
Architecture
CPU - GPU : une comparaison
D'un point de vue simplifié, les processeurs multicores sont de trois types. MIMD (Multiple Instruction, Multiple Data), SIMD ou SPMD (Single Instruction/Program, Multiple Data).Un CPU multicores dispose sur le même socle de quatre cores, de un ou plusieurs caches internes, de caches externes, d'un accès à la RAM par un bus. Les CPU vont par paires dans les Xeon Nehalem pour un total de huit cores. Les CPU sont optimisés pour des accès au mémoires caches à faibles latences, et la logique de contrôle autorise des exécutions spéculative ou dans le désordre.
Là où un processeur générique entre dans le modèle de Von Neumann exécutant une opération à la fois - plusieurs, pas plus de 4 - dans les processeurs superscalaires2.1, les processeurs graphiques entrent dans la classe des ordinateur parallèles : une passe de calcul exécute un programme unique (un kernel ou shader) sur tous les éléments d'un flux de données.
| Processeur | Intel Core 2 Quad | NVIDIA 8800 GTX |
| Paradigme | MIMD | SPMD |
| Cores de calculs | 4 cores | 128 ALUs |
| Fréquence d'horloge | 2.4 GHz | 1.35 GHz |
| Gravure | 65nm | 90nm |
| Millions de transistors | 582 | 681 |
| Cache L1 | 8MB | 16 8K |
| Mémoire | 4GB | 768 MB |
./
En examinant de plus près la distribution des transitors d'un GPU, on peut voir que beaucoup de transistors sont dédiés aux Unités Arithmétiques et Logiquess (ALU), regroupés en multiprocesseur (Streaming Multiprocessors=SM). Chaque SM possède sa propre unité de contrôle et de cache, de dimensions nettement plus réduites que sur des cores de CPU. Les données à l'intérieur du cache d'un SM sont exclusives à celui-ci et ne peuvent être partagées avec les autres. Autrement dit, les flux de données sont indépendants. L'unité de contrôle est moins sophistiquée. En revanche l'intensité de calcul est très haute ; Le flot d'exécution parallèle tire parti d'une bande passante élevée et les latence mémoires entre la RAM embarquée sur le GPU et les caches peuvent être masquées par du calcul, les ALU pouvant travailler pendant que la mémoire est transférée. Dans la pratique, trouver la taille optimale des flots et l'équilibre entre transfert, calcul et restitution nécessiteront souvent le chronométrage des trois opérations.
La carte de calcul Tesla T10
Le carte de calcul Tesla C1060 ou T10 (photo 2.2) est la première carte graphique dédiée au calcul haute performance, conçue pour traiter des problèmes massivement parallèles. Il s'agit d'une version modifiée de la GeForce GTX280, avec davantage de mémoire et pas de sortie vidéo.Processeur hétérogène (photo 2.4), programmable selon un modèle SIMD, une GTX280 possède un socket unique, où sont disposés ses 30 multiprocesseurs (SM) de classe GT200.
La vue ci-dessous schématise l'organisation d'une GTX280. Les SM regroupés par 3 et partagent un niveau de cache.
|
Sa puissance crête atteint environ 400 Gflop/s en précision simple et 78 Gflop/s en double. Le paradigme d'accès parallèle à sa mémoire primaire est multithread avec fusion (ou coalescence) des accès. La bande passante interne est d'approximativement 73,6 Go/s. Le ratio crête flop/octet en précision double est donc proche de 1.
Les bandes passantes externes sur un bus PCI express 16 pour accéder à la RAM du système hôte sont respectivement 8 Go/s et 4 Go/s, pour des zones mémoires alignées et non paginées. Elle chute à 2.5 Go/s en lecture et 2.1 Go/s en écriture sur la RAM du système hôte pour des mémoires non-alignées ou paginées.
La T10 embarque 4 Go de RAM, organisées en huit banques, ces mémoires sont de catégorie GDDR3, cadencées à 1100MHz. Les accès mémoires sont uniformes, à lecture concurrente, écriture aléatoire : plusieurs threads peuvent être ordonnancés pour lire une même zone mémoire, mais que le comportement d'écritures concurrentes est imprévisible.
Selon la charge et les estimations, leur puissance électrique de fonctionnement varie entre 450 (puissance crête) et 235 Watts.
Le multiprocesseur GT200 est un microprocesseur fabriqué avec une précision de 65nm. De type multithread SIMD. Cadencé à environ 1.3GHz, il délivre 2.6 Gflop/s et embarque 64Kb (16 384 registres de 32 bits) de mémoire partagée entre ces 11 unités arithmétiques et logiques : les 8 unités scalaires de précisions flottantes simples (SP), l'unité scalaire de précision flottante double (DP) et les deux unités de fonctions spéciales (SFU).
Chaque unité scalaire exécute jusqu'à 128 threads corésidents partageant un fichier de registre de 2048 entrées. Au total, ce GPU peut exécuter 30720 threads concurrents.
Héritant d'une architecture conçue pour le traitement vidéo et le scintillement exceptionnel d'un pixel isolé étant un phénomène négligeable, les mémoire GDDR2 et 3 des GPU sont dépourvues de code correcteur d'erreur (ECC),
Sur la quantité totale de mémoire physique, environ 50 MB ne peuvent pas être adressés, (la zone mémoire < 0x3A9FFF0), où sont stockés les registres locaux, la table d'allocation des pointeurs et probablement des pilotes.
Le précèdent GPU de Nvidia, GeForce 8800 GTX, dispose de 16 MPs, 128 SP, chacun pouvant gérer 96 threads concurrents, 1024 registres par MP et supporte au maximum 12288 threads concurrents. Conçu à l'échelle de 90nm, il possède 681 millions de transistors et embarque 768 MB de ram DDR2. Chaque processeur de flux dans une GeForce de série 8 peut gérer 96 threads concurrents. Les GT200 sont embarqués sur les cartes GeForce GTX 260, 280, 285, 295, Quadro FX 4800 et 5800 et Tesla C1060.
Fermi/GT300
Les spécifications de la prochaine génération de GPU Nvidia dédiés au HPC sont annoncées en 2010. [7,8,9,10].
Le GPU fermi est composé de 3 milliards de transistors gravés à une échelle de précision de 40nm, répartis en 16 SMP de classe GT300. Comparativement aux Tesla, le nombre de SMP diminue car la quantité d'ALU sur chaque quadruple.
La mémoire globale est de catégorie GDDR5. L'espace d'adressage est en 64-bits pour le système, 40-bits pour l'utilisateur. L'espace mémoire allouable est donc de l'ordre du terabyte, organisé en 6 banques, et cadencés à 3 GHz. La bande passante interne serait de 192 GB/s, sur une interface mémoire de 512-bits. Ces mémoires posséderont des EEC par redondance des transmissions selon le modèle SECDED2.2.
Ce GPU possède un niveau de cache L2, de 768Ko, situé sur la partie centrale du chipset et un système d'adressage et de mise en cache automatique des données de petites tailles depuis l'espace mémoire hôte.
Les mémoires caches embarqués par les SMP passent à 64 Ko, qui pourront être répartis entre le cache de texture et les mémoires partagées sous une configuration 48-14 ou 16-48. Leur fichier de registre est de 32,768 x 32 bit.
Le multiprocesseur GT300 porte 32 unités d'instructions entières et flottantes 32 bits, organisés en deux groupes de 16, plus 16 unités de lecture/écriture (LD/ST) et 4 unités de fonctions spéciales. Ils sont cadencés à 1,5GHz.
L'ordonnanceur global gère 1,536 threads (soit 48 par SMP, correspondants au 2x16 core + les unités LD/ST)
La norme de précision flottante serait améliorée ; Ainsi le résultat intermédiaire des FMA (Fused Multiply Add) aurait une mantisse de 106-bits et 161 bits pendant l'addition.
Les performances annoncées sont respectivement 256 et 512 FMA en précisions simples et doubles par cycle d'horloge, soit 750 Gflops crête en précision double.
Le bus PCIe serait aussi, transformant les modes d'accès bi-directionnel en canal unique, multipliant ainsi par deux la bande passante externe.
La capacité d'exécution concurrente de plusieurs kernels classe ce GPU dans la catégorie MPMD.
De la perspective du programmeur
Il y a quelques années, des programmeurs pionniers découvrirent que les GPU pouvaient être employés pour d'autres traitements que le graphique [11]. Toutefois leur modèle de programmation improvisé était bancale et le pixel shader programmable n'était pas le moteur idéal pour du calcul générique du fait de la liaison avec la partie fixe sur le hardware. Nvidia a évolué ce modèle pour en proposer un meilleur. Le streaming multiprocessor remplace le shader, amalgame qui désignait à la fois le core graphique et le programme correspondant.La programmation des GPU emprunte deux axes : 1) Par programmation directe en utilisant les modèles de programmation délivrés par les constructeurs ou le standard du libre. 2) Par la traduction des programmes existants au moyen de compilateurs appropriés et d'annotations de codes sources par des directives. Le tableau suivant synthétise l'ensemble des solutions retenues pour GPU et processeurs multicores.
| Solution | Spécifique | Multi-OS | Multi-Cores | Compilateur |
| Nvidia CUDA/PTX | NV | ♦ | ||
| OpenCL | ♦ | ♦ | ||
| AMD Brook+/CAL/CTM | ATI | ♦ | ||
| Nvidia Cg + OpenGL | ♦ | |||
| GLSL / OpenGL2+ | ♦ | |||
| Nvidia Cg + DirectX | ||||
| DirectX11 Compute | ||||
| Lib Sh | ♦ | ♦ | ||
| RapidMind | ♦ | ♦ | ♦ | |
| OpenMP | ♦ | ♦ | ♦ | |
| PGI9 | NV | ♦ | ||
| Caps HMPP | ♦ | ♦ |
Les modèles reposants sur des API graphiques ne sont plus considérées comme des solutions viables pour le HPC, mais leur modèle présente un intérêt certain pour les développements de langages de programmation spécialisés ; DirectX11 Compute vise l'accélération des systèmes d'exploitation et des stations de travail ; HMPP et PGI9 sont deux compilateurs concurrents ; OpenMP concerne ici le pilotage des configuration multi-GPU.
Parmis ces solutions, nous retenons pour la suite PGI9, HMPP, CUDA et OpenMP.
Shaders programmables
Avant l'apparition du modèle d'architecture unifiée (programmable par CUDA et OpenCL), le GPGPU a été limité au domaine des API graphiques. Le mode d'adressage était limité aux tailles et aux dimensions des textures, les shaders avaient une capacité de sortie limitée, il n'y avait pas de set d'instruction pour des opérations atomiques entières ou binaires, les échanges d'information restaient limités entre les pixels. L'opération directe de diffusion (scatter) ne peut pas être réalisées avec un processeur de fragment du fait de la localisation fixe de chaque fragment sur la grille lors de sa création. Une opération de diffusion logique implique alors une opération de collecte (gather) additionnelle sur l'intégralité de ensemble.
Toutefois, les langages de shaders tels que GLSL et Cg présentent les avantage d'être simples, compacts, dans un style proche du C++, ils proposent des fonctions intrinsèques de manipulation de matrices et de vecteurs et la composition fonctionnelle. Une série de registres fixes sont préalablement définis, qui correspondent à des zones mémoires segmentées et aux registres qui les accompagnent : tampons de géométrie, tampons de couleurs et d'images, matrices de transformations. Suivant le modèle classique de programmation des GPU avant unification, et conçus pour le rendu d'images, les langages de shaders se composent en règle générale de deux programmes ; un programme de traitement des sommets et un programme de traitement des fragments de pixels. Ceux-ci doivent gérer la liaison avec le pipeline fixe en déclarant quelles sont les variables locales au shader, et quelles sont celles transmisess sur le pipeline.
Le programme final compilé, composé des deux, reçoit en entrée les données qui alimentent les tampons : des géométries structurées par sommets indicés, des matrices de transformations, des coordonnées de textures affectées aux sommets, et des couleurs ou des textures associées. Les deux premières servent à composer la géométrie projective perçue sur un écran, tandis que les suivantes travaillent sur la rasterization de l'image pixélisée. L'ensemble des coordonnées ainsi traité est homogène. Le shader étant séparé du programme hôte et du pipeline fixe, il peut être recompilé entre deux passes de rendu.
La limitation historique de ce modèle tient au fait que la composition finale est placée dans un tampon d'image en écriture exclusive pour être restitué à l'écran uniquement, et dont la relecture est coûteuse. De la même façon, les données placées dans les tampons prédéfinis sont en écriture seule. L'évolution des matériels a permis progressivement de dépasser cette limite en réservant des buffers dans la mémoire du GPU sur lesquels il devient possible de lire et d'écrire.
Les shaders présentent un sujet d'étude très intéressants du fait de la présence d'assembleur, de la liaison fixe avec le matériel et d'un compilateur. Les questions pour mener l'étude à leur sujets sont nombreuses : où réside le programme final ? Quel est ce compilateur et où se trouve-t-il ? Peut-on le reproduire ? L'évolution des GPU permet-elle d'élaborer un nouveau langage dans un style proche de celui des shaders ?
Pour une programmation GPGPU rendue encore plus accessible aux scientifiques, on peut imaginer un langage de programmation simple, disposant d'une interface avec les langages existants, permettant de déclarer les liaisons à l'intérieur d'un programme classique, liaisons attribuées par la suite à un programme écrit dans un paradigme SPMD. C'est ce que font les shaders, mais avec des registres fixes. C'est aussi ce que fait CUDA, mais avec une large mesure de complexité supplémentaire.
CUDA
|
Nvidia CUDA (Compute Unified Device Architecture) est un framework logiciel adapté aux architectures que nous avons présenté en section 2.2. CUDA est une abstraction matérielle : NVidia a toujours masqué l'architecture de ces GPU derrière des interfaces de programmation et le programmeur n'écrit jamais directement "dans le métal". L'abstraction offre deux avantages : 1) un modèle de programmation de haut niveau qui isole des aspects difficile de la programmation des GPU, et 2) l'architecture matérielle peut changer, mais l'API demeure en place, ce qui permettra de conserver les codes préalablement écrits.
Nous n'entrerons pas dans les détails de CUDA et insistons pour renvoyer le lecteur aux travaux de Jimmy Petterson et Ian Wainwright ([12]). Nous recommandons la lecture de leur thèse avant celle les guides de programmation CUDA [13] et [14].
OpenCL
OpenCL est un framework pour écrire des programmes tournant sur des machines hétérogènes : multiprocesseurs, GPU, Cell BE, processeur Fusion et autres. CUDA et OpenCL se ressemblent par bien des aspects. Un retour d'expérience détaillé se trouve aussi dans [12].
CAPS HMPP et PGI9
HMPP (Hybrid Manycore Parallel Programming) de CAPS Entreprise [15] et PGI9 de Portland Group prennent une autre approche, celle du compilateur. La solution existe avant l'avènement des GPU, à l'origine pour faciliter la programmation de machines vectorielles telles que des générations de NEC et de systèmes Cray.
Un programme séquentiel classique est annoté de directives en respectant un placement et une syntaxe qui indiquent où et quand les transferts mémoires surviennent, et quelles fonctions traduire dans un langage SIMD cible avec les données d'entrées-sorties correspondantes.
L'avantage de cette approche est que le code n'est pas modifié mais seulement annoté. Compilé avec un compilateur standard, ces directives sont ignorées au même titre que des commentaires. Compilée avec le compilateur désigné, une portion de code, ou codelet, est généré par le compilateur pour la région à accélérer, produite dans un code source qui peut être ensuite modifiée et dans une librairie partagée appelée à l'exécution. Bien entendu certains algorithmes n'ont pas de propriétés concurrentes parallèles et ceux-ci ne devraient pas être candidats au chargement sur un accélérateur vectoriel. Le choix de la librairie dynamique permet également de porter le code sur différente architecture. HMPP compile des codelets pour les GPU d'ATI/AMD, Nvidia CUDA et les processeurs avec des SSE activées. PGI9 ne le fait que pour les machines supportant les drivers CUDA, mais offre la possibilité d'écrire des Kernels CUDA directement en Fortran, rendant le code plus homogène.
Quelques aspects difficiles de la programmation GPU sont pris en charge par ces compilateurs, tels que les padding mémoires et les recouvrements.
Dans le contexte de code industriel existant, l'emploi de ces directives donne une image plus claire des choix de transferts mémoires à réaliser et permettent d'estimer d'un premier jet si le programme est effectivement accéléré, tout en restant portable sur les différentes architectures. Les kernels peuvent ensuite être modifiés ou réécrits s'il est nécessaire de pousser l'étude de l'optimisation plus loin. Selon la complexité de l'algorithme, le code généré demeure relativement bien inteprétable. Le processus de compilation passe du simple au double, et il peut arriver que parfois de légères modifications du code existant soient nécessaires pour que le compilateur gère correctement l'économie des registres et de mouvements dans les caches.
L'intérêt personnel de maintenir un code historique par l'adjonction de directives pour communiquer avec un générateur de code me semble cependant un exercice assez risqué et à l'intérêt très limité de la perspective du programmeur.
Un exemple de code annoté se trouve en section 5.3.2.
Autres solutions
Les travaux qui suivent ont aboutis à un compilateur qu'il est intéressant de remarquer.
L'accélération d'un solveur de chimie cinétique : une intégration numérique calculatoirement intensive de l'évolution d'espèces chimiques et de leur réactions [16].
[17], le modèle d'advection de traceur qui modélise le transport de constituants atmosphériques sous le forçage des champs des prévisions de vents par le modèle Runge-Kutta du WRF. Le benchmark, employé pour valider le hardware, est un schéma d'advection positif-défini de 5 ordre. L'advection est appelée une fois par itération de Runge-Kutta dans le solveur. Chaque dimension est advectée séparément et ajustée par une approximation en différence-finie.
Les codes des deux projets sont distribués sous leurs versions parallèles (OpenMP) originales et accélérées [18,19].
Les travaux de [20], également sponsorisé par NCAR (National Center for Atmospheric Research) sont axés sur le modèle de résolution de la couche nuageuse sur GPU, par sous mailles en volumes finis. Le modèle historique comporte 8K lignes de Fortran, calcule en précision double. Le problème est petit, mais présente une grande intensité de calcul. Le portage de ce code a suivi la démarche adoptée par le groupe cité au-dessus.
Tout ces travaux menés conjointement par NOAA et le WRF ont aboutis à un compilateur Fortran à C à Cuda (F2C), puis Fortran à Cuda (F2ACC) [21]. Il a été laissé supposé qu'ils soient complétés par l'écriture de kernels Cuda en Fortran compilés via PGI9. Le premier compilateur, F2C, exige des modifications du source existant. Le deuxième non, mais il produit parfois des kernels incorrects. Dans tous les cas, le code généré est passé en revue.
D'autre interfaces semblables à CUDA existent telles que RapidMind [22] et Peakstream [23] (acquis et absorbé par Google en 2008). RapidMind est même allé encore plus loin en adaptant sa plate-forme virtuelle sur Cell, ATI et Nvidia, sans utiliser la plate-forme de ce dernier. Citons enfin GPUocelot, un framework de compilation dynamique pour systèmes hétérogènes, dont un backend CUDA. Licence BSD.
Précision numérique
Introduction
Avant de calculer vite, il faut calculer juste. En l'absence d'un critère d'exactitude,
et c'est le cas des standards à virgule flottante, il est nécessaire de connaître le degré
de précision des calculs.
Le calcul scientifique est très demandeur en précision double de 64 bits, en particulier
lorsque sont appliquées à des techniques itératives sujettes à la propagation d'erreurs d'arrondis
ou des modèles à hautes échelles d'amplitudes.
Le calcul généraliste sur GPU est devenu un exercice très attractif en physique, mais, dans ces processeurs et jusqu'à présent, la vitesse a grignoté sur la précision. Introduire davantage d'erreur de calcul en échange de vitesse nous mène à cette question [24] : Puisque la plupart des calculs en virgule flottante ont de toute façon des erreurs d'arrondis, est-ce important si les opérations arithmétique basiques introduisent un tout petit plus d'erreur d'arrondis que nécessaire ? Si on considère le cas de rendu d'images en temps réel, la réponse est non, le matériel spécialisé calculant dans l'espace discret.
En élargissant le spectre dans le domaine de la simulation, nous essayons de donner un début de réponse pour une classe de matériel et une opération par l'étude des performances et des précisions numériques obtenues par la série des GPU Nvidia T10, en précisions simples et doubles, sur la multiplication de matrices par CUBLAS, l'implémentation CUDA du BLAS.
La régularité et la prédictabilité des accès de données, le requis hautement parallèle et la réutilisation des données en cache, placent cette fonction comme un bon candidat à l'évaluation des performances crêtes des GPU.
Les résultats sont comparés à deux implémentation de référence du BLAS, largement distribuées : ATLAS et MKL. Ils montrent que les précisions sont plus basses, d'un ordre variant entre un et deux, en comparaison avec un calcul sur CPU en simple (32 bits) et double (extension à 80-bits) précisions flottantes.
Produit matriciel parallèle
Soit le produit matriciel C = A × B de dimensions N × N . Le produit calculé parallèlement sans bloc requiert un thread par élément du résultat final. A et B sont chargés N fois depuis la mémoire globale. La réutilisation de donnée de ce calcul est en O(n).
Sur la base d'un calcul par bloc, un bloc de thread de taille b gère une sub-matrice Csub de C de dimensions b×b . A et B ne sont chargés que N/b fois depuis la mémoire globale, ce qui représente l'équivalent en terme de bande passante économisée.
Les blocs résident dans les caches de mémoires partagés.
Le kernel, implémenté dans CUBLAS 2.1, écrit par [25], est une stratégie de mise en cache avec une définition explicite de la dimension des caches optimale pour un GPU GTX280. Les spécifications techniques des Tesla sont fermées par le constructeur, mais les informations données par le driver CUDA indique que le nombre de cores et d'unités de caches est équivalent sur une Tesla T10. Ce ne sera toutefois pas le cas sur les architectures futures.
La dimension de Csub est bornée par le nombre de threads par bloc. En conservant les blocs carrés, et le nombre maximum de threads par blocs étant 512, la taille des blocs est en conséquence 16². Chaque thread porte un identificateur unique sous la combinaison de l'identificateur du bloc et de celui du thread à l'intérieur du bloc, identifiant quelle donnée calculer.
Avec A, B et C respectivement de dimensions (m,k), (k,n) et (m,n). Csub est égale au produit de deux matrices rectangulaires : la sub-matrice A de dimension qui porte les mêmes indices de lignes que Csub et la sub-matrice de B de dimension (n,b) qui a les mêmes indices de colonne que Csub. Ces deux matrices rectangulaires sont divisées en autant de matrices carrées de dimension b nécessaires et Csub est calculée comme la somme des produits de ces matrices carrées. Chacun de ces produits est réalisé en chargeant d'abord les deux matrices correspondantes en cache, chaque thread gérant un élément de la matrice, puis chaque thread calcule un élément du produit. Plus précisément, ce ne sont pas deux sub-matrices qui sont chargés, mais des couple matrice et vecteur variant. Les données du vecteur sont conservées dans les registres. Chaque thread accumule ensuite le résultat de ces produits dans un registre puis le résultat est retourné en mémoire globale. Une première proposition d'une telle répartition de charge apparaît dans [26].
Calculé ainsi en utilisant le maximum de registres et de cache peut être vu comme une stratégie plus générale pour diviser le nombre de relectures. L'idée est applicable au cas de notre super-calculateur à plusieurs processeurs : en divisant de très grande matrices entre plusieurs Tesla connectées par PCIe16, voire à un cluster hybride. Toutes les unités disponibles traitent une matrice m × n. C'est coûteux en mémoire et cela conduit à collecter autant de sub-matrices qu'il y a de charges distribuées. Le calcul peut être asynchrone et recouvert, mais ne passe pas à l'échelle.
|
La solution la plus avancée à l'heure à actuelle pour faire fonctionner de concert une grappe de GPU consiste d'abord à les virtualiser en un seul. CUDA Wrapper [27] est une librairie implémentée en pré-chargement forcé, de telle sorte que les appels à CUDA aux allocations matérielles soient interceptés par elle pour avoir quelques avantages : virtualiser les GPU physiques par un mappage dynamique. Le matériel visible à l'utilisateur est un ensemble consistant de matériel physique et une barrière intrinsèque sur le système partagé prévient les saut accidentels d'un appareil à l'autre. Une affinité NUMA peut être établies entre cores CPU and GPU pour une économie de bande passante.
Spécification de la norme flottante
Selon [13], l'arithmétique binaire flottante des T10 est compatible avec la norme [24]IEEE-754, avec plusieurs déviations.- Les modes arrondis ne sont pas dynamiquement configurables et doivent explicitement définis.
- Pas de mécanisme pour détecter une exception.
- Les opérations se comportent en masquant les exceptions et retournent la réponse masquée telle que définie dans l'IEEE-754 en cas d'exception.
- Pour les mêmes raisons, l'encodage des SNaN est supporté, mais il n'est pas signalé.
- La valeur absolue et la négation ne sont pas compatible avec l'IEEE-754 au regard des NaNs.
Conclusion : pas conforme à l'IEEE-754.
Allocation de mémoire non-paginées sur l'hôte
module cuda
use iso_c_binding
interface
! cudaMallocHost
integer(C_INT) function cudaMallocHost(buffer, size)
bind(C,name="cudaMallocHost")
use iso_c_binding
implicit none
type (C_PTR) :: buffer
integer (C_SIZE_T), value :: size
end function cudaMallocHost
end interface
! alloue dans la RAM hote m*m elements de p octets dans le tableau A
! de memoire alignees sur 64 bits, non paginees pour garantir
! la bande passante maximale entre hote<->gpu
subroutine allocateforgpu( A, m, p, cptr )
! les zones de memoires constante, pile ou tas allant
! dans le gpu requierent la reference d'un pointeur
real, dimension(:,:), pointer :: A
integer :: m
integer :: r
type(C_PTR) :: cptr
integer(C_SIZE_T) :: p
r = cudaMallocHost (cptr, m*m*p)
! bind the C pointer
call c_f_pointer (cptr, A, (/ m, m /))
end subroutine allocateforgpu
|
CUBLAS
Dans sa version 2.1, CUBLAS peut substituer toute autre implémentation du BLAS moyennant quelques modifications mineures et avec cependant plusieurs limitations. Cependant, le comportement par défaut de la librairie est insuffisant : il réalise les allocations, les transferts et les désallocations mémoires sur GPU à chaque appel de fonction d'une routine depuis l'hôte. En général, pour une utilisation efficace, seuls les transferts ou une partie des transferts doivent être effectifs. De plus, la mémoire sur l'hôte doit être alignée et non paginée pour atteindre les performances maximales [14]. Ceci conduit à introduire du code dans un programme existant. Ainsi les transferts de bandes de matrices peuvent survenir plusieurs instructions avant l'appel proprement dit au kernel, et les allocations et désallocations mémoires peuvent avoir lieu une seule fois dans un cycle d'initialisation et de finalisation, et gardées réservées. Ceci est renforcé si le calcul a un comportement redondant. L'autre phénomène remarqué est le temps d'initialisation du GPU, qui a un effet particulièrement notable sur les performances. L'exécution au préalable d'un kernel quelconque y remédie.Si le programme est écrit en Fortran, les appels de requête d'allocation à de la mémoire alignée et non-paginées passe par des liaisons en C pour la manipulation de pointeurs. La comparaison des performances entre une liaison directe de CUBLAS et son usage optimal se voit sur les graphes 4.3 et 4.4 : jusqu'à 100 Gflops d'écarts. La distribution de CUBLAS en version 2.1 est incomplète, dépourvue de toute fonctions complexe de niveau 3. Elles se trouvent en partie dans la version Beta 3.0.
Précision numérique
Pour deviner la précision numérique de CUBLAS sur T10, nous calculons la fonction GEMM C = αAB + βC, de matrices carrées de dimensions m, avec B et C égaux à un entier λ positif. Les colonnes impaires de A sont égales à λ et les colonnes paires de A égales à une valeur flottante ε comprise entre 0 et 1. Le résultat est une matrice C constante.
|
Résultant en C = α(λ²+ελ)k/2 + βλ,
Nous calculons depuis l'espace réel vers l'espace flottant, avec une approximation δ, obtenant le résultat a
R → F
a → a
a = a ( 1 + δ) ~ Ψ
epsilon(double précision) = 2.220446049250313E-16
epsilon(simple précision) = 1.1920929E-07
Résultats
Nous prenons λ et ε=10-n avec n=1..N. Nous avons B=C=2. Comme résultat C=4m+2εm + 4. . La dimension de la matrice varie entre 256 et 12000. La précision est la mesure de la différence relative entre le résultat exact calculé sur le CPU en précision double et le résultat fourni par le GPU.Les graphes 4.7, 4.8, 4.9 montrent les précisions des calculs successifs de l'erreur relative, || C - C' || / C comme une fonction de la taille de la matrice, du résultat C arrondi à l'epsilon machine et des calculs des matrices constantes C' par l'implémentation du BLAS choisie, sujette au cumul des erreurs d'arrondis.
La variabilité et l'annulation des précisions est soit due au mode d'arrondis, à la dénormalisation des nombres ou un calcul exact. Le dernier cas ne survient qu'avec ATLAS pour ε=10-1.
GEMM avec MKL montre une grande stabilité dans la précision. ATLAS n'est exécuté que sur un seul core de CPU et minimse davantage l'erreur relative. Toutefois la fluctuation de la précision est en elle-même un facteur d'imprédictabilité, donc d'imprécisions.
Performances
En 32-bits la fonction SGEMM atteint 370 Gflops de performances crêtes. Ceci survient pour des matrices à partir de carré 12K. A partir de carré 4K, le calcul atteint 300 Gflops. C'est un fait connu en GPGPU : l'efficacité est atteinte avec des matrices suffisamment grandes. En dessous de carré 2K éléments, la latence mémoire à plus d'impact sur le calcul et en dessous d'une certaine dimension, des ALUs sont laissées inactives.
L'augmentation de la taille des matrices à moins d'impact sur l'exécution du kernel que sur le transfert des données. L'utilisation de flots, de copie et d'exécution concurrente permet de diviser les matrices A et C en entrée en bande par lignes et de les transférer en asynchrone, tandis que le kernel calcule une autre bande de données.
Le facteur de gain au regard de la comparaison des deux librairies survient en précision simple pour des matrices supérieures à 1000² , atteignant 2.6 à 7000². En double précision, MKL gagne, et les deux solutions convergent vers une puissance maximale de 67 Gflops.
Remarques
Sur des mémoires dépourvues d'EEC, le calcul répété de grande matrices constantes peut fournir, moyennant l'ajout d'un algorithme de scan parallèle, une détection des occurrences de fautes mémoires et donner une estimation du Temps Moyen Entre Défaillance (TMED). Par exemple, si le TMED pour un core isolé est de 6 mois, sur des unités de 240 cores, il ne serait que de 16 heures et 48 minutes (en établissant une distribution indépendante du temps).Au regard des performances, une image plus précise des mécanismes de synchronisation et de transferts mémoires s'avère nécessaire pour les kernels critiques, à travers l'examen de l'assembleur PTX et l'usage de désassembleurs, tels que Decuda [28]. Celui-ci peut fournir davantage d'information pour l'optimisation des mouvements de caches, comprendre les cycles d'horloge et de manière générale une meilleure compréhension du hardware.
Pour finir, il faut considérer ces architectures multicoeurs sous leur aspect hybride et dans leur ensemble pour obtenir le meilleur entre les deux mondes. Les techniques d'autotuning [29] qui ont déjà inspirées les développement d'ATLAS en sont le meilleur exemple. Sur de vrais problèmes, il est crucial d'obtenir la contribution de toutes les unités de calculs, fonctionnants de concerts lorsque c'est possible, et d'utiliser les plus adaptées aux problèmes à résoudre.
Configuration matérielle du test :
CUBLAS 2.1 sur Nvidia GT200, ATLAS 3.9.15 et MKL 10.0 sur CPU Xeon Nehalem 8 cores X5570 @ 2.93GHz
Cas test
Le calcul selon une méthode aux différence finies, de stencil ou d'interpolations linéaires sur des matrices de l'ordre du teraoctet ne tient pas à l'intérieur de la mémoire d'un ordinateur isolé. Ce traitement, parfaitement parallèle, avec une forte localité des données, nous oblige à utiliser un cluster où les noeuds portent plusieurs coprocesseurs algébriques. Le seul moyen de dépasser les hautes latences de bandes passante entre les unités est de recouvrir les calcul par le transferts asynchrones de fragments de données. Selon un modèle similaire à celui qui se trouve dans une architecture CUDA, avec des blocs de threads calculant pendant que d'autres stockent les données, nous montrons que les délais de synchronisation peuvent être réduits en calculant d'abord les cellules aux interfaces. Nous allons détailler ceci étape par étape.
Méthode aux différences finies
Une méthode d'approximation des dérivées premières, secondes et mixtes utilisant l'expansion des séries de Taylor et les polynômes. Point de départ : les équations de conservations en formes différentielles. Le domaine de la solution est recouvert par une grille. En chaque point de la grille est fait une approximation de l'équation différentielle en remplaçant les dérivées partielles par des approximations des termes aux valeurs nodales des fonctions. Le résultat est une équation algébrique en chaque noeud de la grille. La valeur des variables sur ce noeud et un certain nombre de voisins apparaissent comme des inconnues.(la suite de cette rubrique est publiée uniquement dans la version pdf du document).
Le système d'équation algébriques
Une équation par point de la grille portant la valeur des variables en ce point ainsi que les valeurs des noeuds voisins. Cette équation peut inclure des termes non linéaires. Le processus de résolutions numériques implique une phase de linéarisation. Le système est de la forme :
p est le noeud où l'EDP est approximée et l'index l traverse les noeuds voisins impliqués dans l'approximation en différence finie : p et ses voisins réprésente une "molécule de calcul".
Ce système tridiagonal, creux, est noté matriciellement : A Φ = Q.
A est la matrice carrée des coefficients, Φ le vecteur contenant les valeurs des variables aux noeuds de la grille et le vecteur contenant les termes de la parties droite de l'équation 5.2.
 |
Dans la pratique, on ne réalise pas expressement le produit matrice-vecteur mais une convolution en chaque pixel (ou voxel, en 3d) du domaine discret. Des méthodes plus élaborées, calculant sur des maillages 3d déformables ou non structurés montreront un système similaire, pentadiagonal, non régulier, dans lequel le produit est réalisé. Le système étant creux, une représentation complète de la matrice est prohibitf 5.1.
En 2D, Φ représentera une surface, sous forme d'un tableau linéarisé en mémoire. A est l'opérateur de convolution. Il rassemble et pondère les cellules voisines avec la cellule centrale pour réaliser le calcul 5.1. Le rayon du noyau de convolution est égal à l'ordre de résolution des équations.
L'opération Stencil
En traitement d'image, les filtres de convolutions sont des filtres d'images et le Laplacien est employé comme détécteur sommaire de contour. Ici, l'opérateur stencil évalue les sommes des flux (masses, vélocités, énergies) en chaque point discret de la grille.
Listing 5.1 – L’algorithme de calcul du stencil en deux dimensions de rayon 1
do pp=1,p
do j=2,nlocy+1
do i=2, nlocx+1
ll=i+(j-1)*(nlocx+2)
VV(ll,pp) = s(2)*V(ll-nlocx-2,pp)
+ s(4)*V(ll-1,pp)
+ s(5)*V(ll,pp)
+ s(6)*V(ll+1,pp)
+ s(8)*V(ll+nlocx+2,pp)
end do
end do
end do
|
Plusieurs cas de figure :
Dans le premier cas, on suppose que la charge de calcul est distribuée entre les noeuds. Les interfaces sont échangées à chaque itération sur un canal de communication de haute latence.
Dans le deuxième cas les calculs restent dans le GPU jusqu'à convergence : soit pour une méthode de résolution par itération jacobienne, soit pour la résolution des équations de Laplace.5.2
Le calcul de stencil à 5, 7, 9 ou 27 points a aussi été éxaminé [35] dans le cas d'utilisation d'une machine hybride pour décomposer et calculer spécifiquement le domaine parallèlement sur CPU et GPU.
Transfert des interfaces
Entre chaque itération, les interfaces (ou halos) sont communiquées aux unités de calcul traitant les sous-domaines partageant une interface.
Reprenons la topologie du domaine sur l'architecture de notre machine à calculer. Le nombre d'interface est fonction du niveau de granularité :
soit le produit du diamètre du stencil par la taille des interfaces par le nombre de segmentations. La distribution de la charge est équivalente dans toutes les unités de traitements. Nous avons R, le rayon du stencil de 1 à 5, un problème de dimension 3, et deux subdivisions, m0=4, entre 4 GPU et m1, le nombre de flots à estimer sur la base de la taille du taux d'occupation des GPU et de la dimension du domaine.
Subdivision cartésienne du domaine en sous-domaine
Subdivision des sous-domaine en bandes
Chargement asynchrone des bandes dans chaque GPU
Subdisivion des bandes en blocs
Tant que la solution n'a pas convergé
Mise en cache des blocs
Exécution des kernels
Calcul de la norme résiduelle (Préfixe parallèle)
Si la solution n'a pas convergé
Copie des interfaces dans la mémoire de l'hôte
Section critique : transit des interfaces
Chargement des halos dans les GPU
fin
|
Avec une restriction du problème à un noeud équipé d'une Tesla S1070 (soit 4 Tesla C1060 sur un même bus) le processus d'itération Jacobienne distribué est proposé ainsi :
Ceci est un cas d'école de traitement parallèle. Il est préférable que la dimension des blocs soit le plus uniforme possible. Les délais de communication sont l'étape de synchronisation entre les unités de traitements. La durée du temps de synchronisation est majorée par le temps de transmission de l'interface la plus importante.
Remarque : les changement de contextes pour le pilotage de plusieurs GPU à partir d'un seul noeud sont également en environnement parallèle (pthread ou OpenMP) et la situation de concurrence sur la bande passante lors de l'échange des interfaces est une section critique.
Kernel CUDA
Chaque bloc de threads CUDA traite traite à son tour un sous-ensemble du domaine. La taille des blocs doit être très largement inférieure à la dimension totale de la grille, afin d'avoir suffisemment de blocs pour maintenir tous les multiprocesseurs occupés.
Dans le cas 3D, la quantité limitée de mémoire cache oblige à répéter plusieurs fois le placements de 3 tranches de données du bloc en cache. Le kernel a donc une boucle extérieure pour mise en cache avant synchronisation.
A chaque itération, la tranche inférieure n'est plus nécessaire, la tranche supérieure est chargée, et les pointeurs associés décalés avant le calcul (fig 5.5).
L'optmisation de ce kernel concerne :
- La fusion des accès mémoires de chargements et de stockages
- Le recouvrement entre les warps et/ou blocs actifs pour limiter la latence mémoire
- La minimisation du nombre d'opérations entières d'indexation des tableaux
- L'utilisation maximale du cache et des registres
Ce kernel est encore largement sous-optimal :
Il est prévu pour des domaines de toute dimension, et la grille étant dynamiquement configurée, le module peut ne pas être entier : les blocs de threads doivent déterminer s'il se situent sur la partie non-entière d'un module et à l'intérieur de ce bloc, tester s'ils sont encore à l'intérieur du domaine. Or comme on le sait, les GPU n'aiment pas les branchements conditionnels.
Les accès des halos latéraux sont non fusionnés, divergents, les conflits de banques multiplient par 2x les lectures mémoires.
Les 16 Ko de mémoire cache sont le facteur limitant à la dimension des blocs : dépasser cette limite impose une stratégie de mise en cache circulaire. Cependant, l'absence de contrôle de l'ordonnanceur, excepté la barrière de synchronisation, rend la performance du globale kernel relative au thread le plus long. Il convient donc de réduire le nombre d'itérations individuelles des threads 5.3.
|
L'histogramme 5.6 souligne le goulet que représente la bande passante et l'importance de minimiser ces transferts.
Relaxation des synchronisations
Reprenons le cas d'utilisation d'une ou de plusieurs Tesla S1060 (4 GPU). Une relaxation des synchronisations peut être envisagée pour que l'échange des donnés recouvre les calculs. On peut également mettre en oeuvre une barrière de synchronisation fine pour que, moyennant un coût mémoire supplémentaire, plusieurs itérations successives soient réalisés sur un domaine auquel est soustrait un halo interne à chaque itération. Pour chaque itération, le halo soustrait qui est généré autour du domaine restreint correspond au sous-ensemble non calculé et en attente des halos voisins. Ceci introduit une nouvelle difficulté sous la forme d'une subdivision moins uniforme : la texture représentant le domaine prend la forme d'une pyramide discrète, dont il faut calculer le complément. Cette solution maximize l'utilisation du cache et occupe plus d'espace mémoire. Les itérations successives respectent une compléxité logarithmiques. En revanche le calcul des compléments introduit davantage de cellules fantômes.
Calculer les cellules fantômes en premier
Revoyez la figure 5.4. A une échelle gros grain, la barrière de synchronisation entre les noeuds ne tient pas au calcul de l'intégralité des sous-parties, mais uniquement des interfaces. Il faut donc commencer par calculer d'abord les interfaces. Ceci fait, procéder immédiatement à la synchronisation.
Sur l'accélération des codes
A partir du niveau 3, et en simple précisions, les fonctions BLAS sont plus rapides sur GPU à condition que les dimensions de matrices soient suffisantes. A cet effet, nos benchmarks sur les précision numériques (4.1) peuvent fournir les rapports de vitesse existants entre CUBLAS et MKL selon les dimensions des matrices et la précision pour toutes les fonctions du BLAS 2 et 3.Cet exercice présente un intérêt à la fois technique, académique et bref.
Une librairie, GBLAS, a été instruite pour fournir les interfaces C d'appel aux fonction BLAS que peuvent invoquer les programmes écrits en Fortran. Selon les scénarios résultants des benchmarks évoqués ci-dessus, le code final devra faire appel aux fonctions présentes dans l'une ou l'autre des librairies.
In fine, une librairie hybride, destinées à supplanter LAPACK, avec des capacités d'auto-calibrage en compilation et en exécution pour pouvoir déterminer le choix précédent. C'est la technique mise en oeuvre dans ATLAS. L'objectif est de garder cette librairie adaptable aux changements rapides de matériels. Le projet MAGMA s'y emploie, développé à l'ICL [36,37,38]. L'évolution matérielle modifie le comportement optimal d'un kernel car celui-ci est lié aux quantités de mémoires partagées et de registres disponibles [29].
Le calcul des valeurs propres est un problème qui se prète bien au traitement sur unités vectorielles [39],[40].
Les résultats en précision double ont montré que les performances crêtes obtenues sur une Tesla était au dessous des bi-quad cores. Ceci ne tient qu'au nombre de FPU des Tesla, et le rapport 1/8 unités SP et DP devrait passer à 1/2 dans les architectures Fermi.
Axes de recherches sur la représentation des maillages non-structurés en carte graphique :
- Commencer par calculer les normales des faces.
- Rechercher une représentation compacte en mémoire qui réduisent les caches-manqués et puissent tenir dans les caches des multiprocesseurs.
- Algorithmes de coloration de graphes [41].
Pour finir, des travaux sur la transformée rapide de Fourrier obtiennent des gain d'accélération remarquables (jusqu'à 150x pour un cas de détection radar [12]), mais avec des précisions numériques très légèrement inférieure (ce qui importe peu dans les traitements de signaux en temps rééel d'un signal d'office bruité).
Conclusions
Les nouveaux modèles de programmation des GPU n'ont plus grand chose à avoir avec la programmation GPGPU via des API graphiques.
Nous avons évoqué l'idée d'un langage qui aurait l'avantage de la simplicité des shaders. Des kernels écrits en Fortran pourraient bien remplir ce rôle.
Si nous avons comparé GPU et CPU, on peut aussi constater que l'architecture des cartes graphiques unifiées s'assimile peu à peu à un micro-cluster. La principale différence depuis la perspective du programmeur devient le contrôle de la mémoire très explicite, et l'adaptation des problèmes de calculs aux contraintes de ce hardware simple, mais puissant. La plupart des difficultés avec la programmation parallèle, telle que le passage des messages et les deadlocks, n'existent pas dans CUDA pour la simple raison que le passage de message est impossible à la façon des grappes de CPU. Un mécanisme de synchronisation globale et locale le remplace.
Le bénéfice, au delà de plus de puissance de calcul est que les problèmes associés aux fonctionnalités mentionnées manquantes n'arrivent jamais : l'architecture interdit tout simplement le programmeur de tomber dans un tel piège.
Les cartes GPU sont peut-être amenée à disparaître, comme disparurent les cartes sons, et à se retrouver intégrés dans les CPUs. Le modèle de programmation restera. A terme, on peut supposer que les compilateurs et les processeurs seront devenus suffisamment sophistiqués pour réinterpéter automatiquement la logique des codes historiques en termes d'exécution parallèle.
Annexe
Liste des GPU compatibles CUDA
et capacité de calcul 1.3
Les cartes compatibles avec CUDA et de capacités de calcul 1.3
implémentent l'ensemble des spécifications de CUDA et peuvent calculer en double précision IEEE-754 déviée.
Elles sont toutes basées sur le chipset GT200, exceptée les chipset T10 de Tesla, qui sont une version modifiée des GT200.
le multiplicateur indique le nombre de GPU embarqué dans la carte.
 |
Tesla S1070
ne S1070 est un {\em rack} de 4 TeslaC1060. Le GPU le plus puissant avant l'arrivée du Fermi.
Présente un intérêt pour le HPC qui nécessite impérativement une solution redondante ou les 4Go de Mémoire qu'elle embarque.
|
 |
Tesla C1060
Une solution HPC à l'intérieur d'une station de travail. Deux cartes autorise à débuter en programmation environnement Multi-GPU, ce qui présente un challenge.
Consommation électrique : 150 W. Fréquence des SPU : 1.296 à 1.44 Ghz
|
 |
Quadro Plex 2200 D2
Sous forme de station de travail ou de rack. Cher. 4 port DVI. jusqu'à 8 écran peuvent s'y brancher. Spécialisé pour le rendu d'image complexe ou interactif.
|
 |
Quadro FX 5800
Un GPU pour le calcul ou la visualisation. A comparer à une Tesla mais avec une sortie vidéo. Ce qui facilite et rend plus agréable la programmation.
|
 |
Quadro FX 4800
Similaire à une 5800 mais avec des capacités moindres.
|
 |
GeForce GTX 295
Deux GTX200 à l'intérieur d'une seul carte. Le système verra deux cartes.
Nécessite une alimentation électrique d'au minimum 680W. Recommandé pour
les phases de prototypages.
|
 |
GeForce GTX 285
Alim. Électrique 550 W minimum. Fortement recommandé pour la phase de prototypage.
|
 |
GeForce GTX 280
Fréquence d'holroge processeur (MHz) 1242 MHz
|
 |
GeForce GTX 260
Presque identique à la 280.
|
Bibliography
- 1
-
André Lebeau.
L'engrenage de la technique.
Gaillmard, 2005. - 2
-
Leo Gerat.
Le calcul scientifique se met au vert, Juillet-Août 2009. - 3
-
Herb Sutter.
The Free Lunch Is Over: A Fundamental Turn Toward Concurrency in Software, 2005. - 4
-
Jon Peddie.
GPU market defies gravity so far, April 2009. - 5
-
CUDA Zone.
http://www.nvidia.com/cudazone. - 6
-
Michael Garland, Scott Le Grand, John Nickolls, Joshua Anderson, Jim Hardwick,
Scott Morton, Everett Phillips, Yao Zhang, and Vasily Volkov.
Parallel Computing Experiences with CUDA.
IEEE Micro 28, August 2008. - 7
-
Tom Halfhill.
Looking Beyond Graphics.
Technical report, In-Stat, 2009. - 8
-
David Patterson.
The Top 10 Innovations in the New NVIDIA Fermi Architecture and the Top 3 Next Challenges .
Technical report, Parallel Computing Research Laboratory, U.C. Berkeley1, September 2009. - 9
-
Glaskowsky.
NVIDIA’s Fermi: The First Complete GPU Computing Architecture.
Technical report, In-Stat, 2009. - 10
-
Nathan Brookwood.
NVIDIA Solves the GPU Computing Puzzle.
Technical report, September 2009. - 11
-
Gpgpu.org.
http://www.gpgpu.org. - 12
-
Jimmy Pettersson and Ian Wainwrigth.
Radar Signal Processing with Graphics Processors.
January 2010. - 13
-
NVIDIA CUDA - Programming Guide.
Technical report, NVIDIA, 4 2009.
Version 2.3. - 14
-
NVIDIA CUDA - Programming Best Practices Guide.
Technical report, NVIDIA, July 2009.
Toolkit v2.3. - 15
-
Romain Dolbeau, Stéphane Bihan, and François Bodin.
HMPP : A Hybrid Multi-core Parallel Programming Environment.
October 2007.
Workshop on General Processing Using GPU, Boston, October 2007. - 16
-
J. Linford, J. Michalakes, A. Sandu, and M. Vachharajani.
Multi-core acceleration of chemical kinetics for simulation and prediction.
2009. - 17
-
John Michalakes and Manish Vachharajani.
GPU Acceleration of Numerical Weather Prediction.
2009. - 18
-
John Michalakes and Manish Vachharajani.
GPU Acceleration of Scalar Advection.
http://www.mmm.ucar.edu/WG2bench/. - 19
-
John Michalakes, John Linford, , Manish Vachharajani, and Adrian Sandu.
GPU Acceleration of a Chemistry Kinetics Solver.
http://www.mmm.ucar.edu/wrf/WG2/GPU/Chem_benchmark/. - 20
-
Evan F. Bollig.
Gpu acceleration of the cloud resolving.
Technical report, 2009. - 21
-
UCAR.
F2C.
http://www-ad.fsl.noaa.gov/ac/Accelerators.html. - 22
-
Writing Applications for the GPU Using the RapidMind Development Platform.
Technical report, Rapidmind, 2006. - 23
-
Matthew Papakipos.
The PeakStream Platform : High-Productivity Software Development for multi-Core Processors.
Technical report, April 2007. - 24
-
David Goldberg.
What every computer scientist should know about floating-point arithmetic, volume 23.
ACM, New York, NY, USA, 1991. - 25
-
Vasily. Volkov and James Demmel.
LU, QR and Cholesky factorizations using vector capabilities of GPU.
Technical Report No. UCB/EECS-2008-49, June 2008. - 26
-
K. Fatahalian, J. Sugerman, and P. Hanrahan.
Understanding the Efficency of GPU Algorithms for Matrix-Matrix Multiplication.
Graphics Hardware, october 2004. - 27
-
National Center for Supercomputing Applications University of Illinois,
Innovative Systems Lab.
Cuda Wrapper.
http://sourceforge.net/projects/cudawrapper/. - 28
-
Wladimir J. Van Der Laan.
Decuda.
http://wiki.github.com/laanwj/decuda. - 29
-
Yinan Li, Jack Dongarra, and Stanimire Tomov.
A Note on Auto-tuning GEMM for GPU.
January 2009. - 30
-
Jeff Bolz, Ian Farmer, Eitan Grinspun, and Peter Schröoder.
Sparse matrix solvers on the GPU: conjugate gradients and multigrid.
In SIGGRAPH '03: ACM SIGGRAPH 2003 Papers, pages 917-924, New York, NY, USA, 2003. ACM. - 31
-
Luc Buatois.
Algorithmes sur GPU de visualisation et de calcul pour des maillages non-structurés.
PhD thesis, Institut National Polytechnique de Lorraine, 2008. - 32
-
Luc Buatois, Guillaume Caumon, and Bruno Lévy.
Concurrent Number Cruncher, A GPU implementation of a general sparse linear solver.
2007. - 33
-
Nathan Bell and Michael Garland.
Efficient Sparse Matrix-Vector Multiplication on CUDA.
December 2008. - 34
-
Eric Darrigrand, Sébastien Giraud, and Vincent Pit.
Une méthode multipôle rapide pour l’équation de Laplace 3D. - 35
-
Sundaresan Venkatasubramanian.
Tuned and asynchronous stencil kernels for CPU/GPU systems, May 2009. - 36
-
Vasily Volkov and James Demmel.
Benchmarking GPU to tune dense linear algebra.
November 2008. - 37
-
Marc Baboulin, Jack Dongarra, and Stanimire Tomov.
Some Issues in Dense Linear Algebra for Multicore and Special Purpose Architectures.
2008. - 38
-
Marc Baboulin, Jack Dongarra, and Stanimire Tomov.
Towards dense linear algebra for hybrid GPU accelerated manycore systems.
2008. - 39
-
Vasily Volkov and James Demmel.
Using GPU to Accelerate the Bisection Algorithm for Finding Eigenvalues of Symmetric Tridiagonal Matrices.
Technical Report UCB/EECS-2007-179, EECS Department, University of California, Berkeley, Dec 2007. - 40
-
Chritian Lessig.
Eigenvalues computation with CUDA, October 2007. - 41
-
K.A. Hawick, A. Leist, and D.P. Playne.
Parallel Graph Component Labelling with GPU and CUDA.
June 2009.
About this document ...
Footnotes
- ... superscalaires2.1
- Les accès par cache du Pentium 4 lui permettent de réunir une valeur SSE de 128 bits (un paquet de 4 flottants de 32 bits) en un cycle.
- ... SECDED2.2
- single (bit) error correction, double error detection
- ... prohibitf5.1
- Ceci est renforcé s'il faut faire des transferts mémoire depuis l'hôte vers les unités de calcul. Le produit vecteur - matrice creuse fait l'objet de nombreux travaux dans le monde du calcul scientifique. Le problème est bien adapté au GPU : traité sur des architectures non-unifiés avec des méthodes de GPGPU traditionnelle [30] ou sur des architectures plus récentes (CUDA et CTM) dans [31,32] pour la méthode du gradient conjugué préconditionné par Jacobi. Plusieurs formats on été testés sur des variétés de matrices courantes dans leur domaine d'application[33]. Les meilleurs cas atteignent respectivement 15 et 10 Gflops, en précisions simple et double.
- ... Laplace.5.2
- Une résolution spécifique de l'équation de Laplace en 3D éxige des méthodes numériques plus musclées pour réduire la compléxité en temps et en mémoire. Multipôle rapide [34] et convolutions/déconvolutions par FFT.
- ... threads5.3
- Nota : pas dans un cas général sans branchement conditionnel et où tous les accès seraient alignés.
