Introduction à la programmation GPU - Part 3
Introduction à l’architecture CUDA
Dans la première partie de cette série consacrée à l’introduction à la programmation GPU (Graphics Processing Unit), nous avons présenté les différentes technologies et l’algorithme qui sera utilisé pour chaque technologie. Nous avons constaté les performances de cet algorithme à la fois en mode séquentiel et en mode parallèle sur un CPU (Central Processing Unit) de type quad-cœurs. Dans la seconde partie nous avons présenté comment installer l’environnement de développement CUDA C sous Windows avec Visual Studio 2010. Nous avons aussi montré comment interroger une carte graphique compatible CUDA afin de récupérer ses caractéristiques techniques.
Dans cette troisième partie, nous allons vous présenter les fondamentaux de l’architecture CUDA afin de vous sensibiliser à la programmation GPU qui réclame une attention importante sur l’architecture matérielle. Nous allons pour cela implémenter avec CUDA l’algorithme présenté dans la première partie.
Architecture CUDA
Si vous êtes néophyte en programmation CUDA, comprendre l’organisation des threads GPU d’une architecture CUDA est essentiel pour commencer à programmer vos propres algorithmes. Vous devez considérer les threads GPU comme à la fois plus nombreux, plus légers et plus rapides à lancer que leurs homologues CPU. Par exemple, la contention mémoire que peut provoquer un nombre très important de threads CPU a beaucoup moins de sens avec des threads GPU. Les threads GPU sont organisés d’une manière totalement différente de ceux d’un CPU. Les threads GPU sont répartis dans des structures de type blocs, qui sont elles-mêmes étagées dans une structure de type grille, figure 1. De la même manière, plusieurs grilles peuvent coexister si vous possédez plusieurs cartes graphiques.
Figure 1- Schéma simplifié de l’architecture physique GPU
Les blocs de threads représentés ici par les lignes horizontales contiennent sept threads chacun. Il y a six lignes de blocs de threads sur la figure 1. Sur le plan de la mémoire, chaque bloc de threads possède un contrôleur de registres (rectangle orange) et une zone de mémoire partagée (rectangle rouge) accessible par tous les threads du bloc. Cette mémoire partagée est extrêmement rapide et elle peut être utilisée pour accélérer les performances de vos algorithmes. Une grille contient aussi de la mémoire globale (rectangle rouge DRAM) plus lente que la mémoire partagée, très simple à utiliser, mais surtout bien plus volumineuse. C’est à travers cette mémoire que les échanges de données se font entre votre host (thread CPU responsable du lancement du traitement GPU) et votre carte graphique. Le host alloue la mémoire nécessaire sur le GPU, puis transfère les octets recopiés dans la mémoire allouée précédemment, mémoire host vers mémoire GPU et inversement dans le cas d’un retour du code GPU. Nous reviendrons plus en détail sur l’organisation mémoire dans un prochain article, mais pour l’instant passons à l’organisation des threads GPU.
Organisation des threads
Généralement, l’exécution d’un programme CUDA, se traduit par le lancement d’une ou plusieurs fonctions dites « noyau ». Ce nom peut vous sembler très technique, mais il signifie simplement que le code respectif de ces fonctions va s’exécuter sur la carte graphique et non sur les cœurs CPU de votre machine. Ces fonctions sont en charge de définir le nombre de threads souhaités pour l’exécution. Dès qu’une fonction noyau est lancée, une organisation de threads se met en place pour satisfaire la demande de calcul. Cette organisation se traduit par une grille contenant souvent des milliers voire des millions de threads.
Figure 2 - Organisation des threads GPU sur une architecture CUDA
Cette organisation physique peut vous sembler étrange, mais vous devez vous familiariser avec l’organisation atypique des threads GPU que l’on retrouvera sous une forme quasi identique en OpenCL ou en C++ AMP. En effet, pour programmer votre calcul GPU, vous allez devoir vous familiariser avec ce découpage structurel qui vous permettra progressivement de vous repérer facilement. Vous êtes en charge d’exprimer combien de blocs et threads sont nécessaires pour exécuter votre algorithme. Sur la figure 2, la grille verte contient six blocs indicés via un numéro de ligne et un numéro de colonne. C’est ainsi que nous repèrerons les blocs dans nos fonctions noyaux. Chaque bloc contient 12 threads, les indices des threads permettent à votre code de savoir où il se trouve vis-à-vis de cette topologie de 72 threads. Enfin, notons que nous aurions pu représenter les blocs sous la forme d’un cube, car le système de coordonnées de blocs repose en réalité sur 3 dimensions.
Si vous souhaitez paralléliser vos codes via CUDA C, vous devez tenir compte de cette organisation. En effet, calculer le nombre threads par bloc et le nombre de blocs est nécessaire pour paralléliser votre traitement. Généralement votre traitement réclamera un nombre d’itérations supérieur à un bloc de threads. Dans ce cas, vous aurez besoin d’utiliser plusieurs blocs pour satisfaire l’exécution de votre algorithme. Dans notre exemple, le produit matriciel de référence engage des matrices de taille 1024 * 1024, ce qui réclame 1 048 576 threads !
L’exécution de l’outil deviceQuery.exe, figure 3, présenté dans l’article précédent nous permet d’obtenir de nombreuses informations sur l’architecture CUDA courante.
Figure 3 - Exécution de l’outil DeviceQuery.exe
Mais dans le cadre de cette introduction, nous allons nous attarder sur les éléments cerclés en rouge sur la figure 3. En effet, lorsque nous lancerons des calculs sur la carte graphique nous serons heureux de connaître le nombre maximum de threads par blocs, la taille maximum d’un bloc et la taille maximum d’une grille. Tout développeur CUDA doit tenir compte de ces informations pour paralléliser ses traitements.
Communication entre le code CPU et le code GPU
Lorsque vous lancez un programme qui réalise un calcul sur GPU, vous utilisez malgré vous un thread système par défaut pour démarrer votre application. Un thread CPU est à la fois nécessaire pour exécuter votre programme, mais aussi pour héberger le code noyau CUDA, avant son transfert sur la carte graphique. Ce code noyau sera exécuté par de nombreux threads GPU, qui participeront massivement à la réduction temps de calcul. Lorsque nous parlons de code GPU, nous parlons d’un code qui s’exécute sur l’architecture CUDA. En effet la programmation GPU distingue deux modes d’exécution : host et kernel. Le mode host, est le mode d’exécution par défaut de tout programme Win32, mais la programmation GPU permet d’exécuter votre algorithme sur votre carte graphique. Pour clarifier l’exécution d’un programme CUDA, je vous propose le schéma de la figure 4.
Figure 4 - Protocole de communication pour un appel noyau, entre le mode CPU/Host et le mode GPU/Noyau
C’est la méthode MatrixMultiplyCudaSquare qui illustre l’essentielle des étapes décrite sur la figure 4. Cette fonction appartient au fichier « MultMatrix.cu » à télécharger avec cet article.
void MatrixMultiplyCudaSquare(float* A, float* B, float* C, int size)
{
float* deviceA;
float* deviceB;
float* deviceC;
size_t sizeAlloc = size*size*sizeof(float);
checkCUDA(cudaMalloc((void**)&deviceA, sizeAlloc), "cudaMalloc");
checkCUDA(cudaMalloc((void**)&deviceB, sizeAlloc), "cudaMalloc");
checkCUDA(cudaMalloc((void**)&deviceC, sizeAlloc), "cudaMalloc");
checkCUDA(cudaMemcpy(deviceA, A, sizeAlloc,cudaMemcpyHostToDevice), "cudaMemcpy");
checkCUDA(cudaMemcpy(deviceB, B, sizeAlloc, cudaMemcpyHostToDevice), "cudaMemcpy");
checkCUDA(cudaMemcpy(deviceC, C, sizeAlloc, cudaMemcpyHostToDevice), "cudaMemcpy");
dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 numBlocks(adjust(size) / threadsPerBlock.x, adjust(size) /threadsPerBlock.y);
MatrixMultiplySquare<<<numBlocks, threadsPerBlock>>>(deviceA, deviceB, deviceC, size);
checkCUDA("MatrixMultiply");
checkCUDA(cudaMemcpy(C, deviceC, sizeAlloc, cudaMemcpyDeviceToHost), "cudaMemcpy");
checkCUDA(cudaFree(deviceA), "cudaFree");
checkCUDA(cudaFree(deviceB), "cudaFree");
checkCUDA(cudaFree(deviceC), "cudaFree");
}
La fonction MatrixMultiplyCudaSquare en charge de lancer l’algorithme noyau
Reprenons chacune des étapes avec une illustration du fragment du code C associé.
Étape 1 - Définir et allouer les paramètres qui devront être copiés côté GPU pour les fournir à la fonction noyau
La transmission des paramètres à notre fonction noyau, doit passer par des variables dont l’allocation mémoire est issue de la carte graphique. Dans notre cas, nous devons recopier les 3 vecteurs représentants nos matrices A, B et C dans des variables préfixées par le mot device, deviceA, deviceB et deviceC. Pour allouer de la mémoire dans la mémoire globale de la carte graphique, nous utilisons la fonction cudaMalloc(), qui reprend le principe de la célèbre fonction malloc() du langage C. Pour calculer la taille de nos matrices (nbCol*nbRow*sizeof(type)), nous multiplions le nombre de colonnes par le nombre de lignes (ici c’est une matrice carrée) et par la taille de l’item, ici un nombre flottant.
void MatrixMultiplyHost(float* A, float* B, float* C, int size)
{
float* deviceA;
float* deviceB;
float* deviceC;
size_t sizeAlloc = size*size*sizeof(float);
cudaMalloc((void**)&deviceA, sizeAlloc);
cudaMalloc((void**)&deviceB, sizeAlloc);
cudaMalloc((void**)&deviceC, sizeAlloc);
Les variables noyau et leurs allocations
On passe à la fonction cudaMalloc l’adresse du pointeur mémoire devant être alloué et la taille du buffer à allouer. Les matrices A et B sont les valeurs d’entrées. La matrice C représente la matrice résultat en sortie.
Étape 2 - Copier les paramètres vers la mémoire GPU
Après avoir alloué les variables GPU, nous devons recopier respectivement le contenu des matrices A, B, C dans la variable deviceA, deviceB, deviceC. Nous utilisons la fonction cudaMemcpy et le drapeau cudaMemcpyHostToDevice afin d’assurer le transfert des données dans le sens CPU vers GPU.
cudaMemcpy(deviceA, A, sizeAlloc, cudaMemcpyHostToDevice);
cudaMemcpy(deviceB, B, sizeAlloc, cudaMemcpyHostToDevice);
cudaMemcpy(deviceC, C, sizeAlloc, cudaMemcpyHostToDevice);
Copier les variables hosts dans les variables noyaux
Arrivé à ce stade, les trois matrices « noyaux » deviceA, deviceB et deviceC sont allouées et initialisées et recopiées dans la mémoire globale GPU.
Étape 3 - Appel de la fonction noyau
Nos variables sont prêtes et transmises au code noyau, nous pouvons donc appeler la fonction noyau. Cependant, nous devons préciser combien de fois la carte graphique devra exécuter notre algorithme. En effet en programmation CUDA C, nous devons préciser combien de fois le code noyau devra être exécuté. Nous touchons un point très important de la programmation parallèle sur GPU. En programmation CUDA C, les boucles sont prises en charge par le moteur d’exécution CUDA et c’est donc le code appelant, côté host, qui se doit de préciser le nombre d’itérations à réaliser.
MatrixMultiply<<<numBlocks, threadsPerBlock>>>(deviceA, deviceB, deviceC, size);
Remarque : La programmation parallèle GPU permet de résoudre des problèmes de performance sur des algorithmes orientés données. D'un point vue de la programmation, cela revient à optimiser les boucles « for » des implémentations et dans notre cas nous sommes dans le cadre de deux boucles imbriquées.
Le lancement de notre méthode noyau doit respecter la syntaxe suivante ; nom de la méthode, suivie de trois chevrons ouvrant, le nombre de blocs nécessaires, une virgule, le nombre de threads par blocs, suivit de trois chevrons fermants. Le reste de la signature ne présente pas de différence avec une fonction C standard. Si vous souhaitez lancer 512 fois votre fonction, vous pouvez par exemple vous exprimer ainsi:
MatrixMultiply <<<1, 512>>>(deviceA, deviceB, deviceC, size);
Remarque : La syntaxe pour appeler une fonction noyau ne peut être reconnue par un compilateur C standard car la syntaxe avec les chevrons n'est pas reconnue. Pour que votre compilateur CUDA reconnaisse ce code, vous devez le placer dans un fichier de type « .cu ».
Dans notre exemple nous souhaitons calculer le produit matriciel de matrices 1024 * 1024. Comment faire pour exécuter 1 048 576 threads exigés par notre traitement. Rappelez-vous de l’organisation des threads sur l’architecture CUDA, nous disposons de nombreux blocs. Dans les faits, nous disposons d’un cube de blocs (1024, 1024, 64) par grille. Pour calculer le nombre d’itérations des deux boucles externes de notre algorithme, nous utiliserons deux variables CUDA de type dim3 supportant une configuration cubique (x, y, z). Dans notre cas nous utiliserons que deux dimensions des structures dim3: threadsPerBlock pour définir combien de threads seront créés par bloc et numBlocks pour définir combien de blocs seront utilisés. Cette technique de répartition entre les blocs et le nombre de threads par bloc permet d’exploiter la topologie de threads CUDA efficacement. Par exemple, nous définissions pour chaque bloc de threads la taille suivante : BLOCK_SIZE * BLOCK_SIZE.
dim3 threadsPerBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 numBlocks(adjust(size) / threadsPerBlock.x, adjust(size)/threadsPerBlock.y);
Déclaration des variables spécifiant le nombre de threads à lancer
Sachant que la valeur de BLOCK_SIZE est de 32 dans notre programme, ce qui signifie que nous définissons un plan de 1024 threads par bloc. Côté blocs, nous définissons un plan de 1024/32 * 1024/32 : 1024 blocs. Finalement, nous obtenons les 1048576 threads attendus. En procédant ainsi, nous répartissons notre demande de calcul sur l’architecture CUDA sans difficulté.
Cependant, l’utilisation de BLOCK_SIZE, peut poser un problème. Pour éviter des soucis, si la valeur size n’est pas modulo de la valeur de BLOCK_SIZE, nous passons par l’intermédiaire de la fonction utilitaire adjust() , qui nous assure un résultat positif, mais parfois un peu supérieur à notre demande.
int adjust(int n)
{
if (n < BLOCK_SIZE) { return n; }
return (n / BLOCK_SIZE + (n % BLOCK_SIZE == 0 ? 0 : 1)) * BLOCK_SIZE;
}
Ajuster le nombre de threads vis-à-vis de la valeur BLOCK_SIZE
Nous verrons, lorsque nous étudierons le code noyau, comment rectifier ce petit excès potentiel dans le nombre d’itérations.
MatrixMultiply<<<numBlocks, threadsPerBlock>>>(deviceA, deviceB, deviceC, size);
Lancer le traitement sur la carte graphique
Remarque : si vous lancez un nombre de threads supérieurs à la capacité de votre carte graphique, le moteur d’exécution CUDA retournera une erreur de type : « invalid configuration argument ». Il est donc important d’être vigilant sur le paramétrage de fonction noyau afin d’éviter de mauvaises surprises.
Après le lancement de la fonction noyau, le code noyau va s’exécuter sur l’architecture CUDA. Noter qu’il n’est plus possible de modifier le nombre de threads une fois, le traitement lancé.
Étape 4 - Exécuter l’algorithme sur l’architecture CUDA de votre carte graphique
Du côté de l’implémentation de notre calcul matriciel, vous allez sans doute noter que le code de l’algorithme a un peu changé.
__global__
void MatrixMultiply(float* A, float* B, float* C, int size)
{
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
if (j < size && i < size)
{
float sum = 0;
for (int k = 0; k < size; ++k) {
sum += A[i * size + k] * B[k * size + j];
}
C[i * size + j] = sum;
}
}
Code de notre algorithme décliné en une fonction kernel
Le mot clef __global__ en préfixe de la définition de la fonction MatrixMultiply « Noyau » permet de signaler au compilateur CUDA, nvcc, que cette fonction est globale et donc appelable depuis le host, mais s’exécute exclusivement sur la carte graphique.
Les extensions du compilateur CUDA C offrent plusieurs types de décorations pour les fonctions. Le tableau ci-dessous, résume les usages possibles des différents types de décoration.
Si votre fonction kernel doit être décomposée en sous fonctions devant être appelable depuis le périphérique, les sous fonctions doivent être décorés avec le préfixe le mot clef __device__ . Pour les fonctions s’exécutant exclusivement sur le CPU, le préfixe __host__ est utilisable, mais pas utile pour compiler correctement.
Contrairement à la version séquentielle où nous bouclions sur les lignes et les colonnes des matrices A et B, l’exécution GPU assure ces boucles en répartissant le code de notre algorithme (code noyau). Sur chacun des threads GPU réclamés (1048 576 threads), chaque instance de la fonction noyau se trouve dupliquée sur chacun des threads impliqués dans notre traitement. Ainsi, chaque code noyau obtiendra automatiquement les informations sur les indices du thread GPU qui l’exécute.
int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;
Révision de nos variables vis-à-vis du moteur d’exécution CUDA
Les variables multidimensionnelles threadIdx, blockDim et blockIdx permettent de retrouver les indices du thread courant, afin de renseigner nos variables d’origines, i et j. Dans notre exemple de produit matriciel, nous disposions de deux boucles externes, indicées respectivement par i et j. On note que suite au calcul sur le nombre de threads réclamé par le code appelant, le code noyau se protège via l’instruction if, pour ne traiter que les indices valides vis-à-vis de la taille de nos matrices.
L’architecture CUDA, permet de répartir le calcul sur une grille de threads GPU. Sur la figure 5, nous comparons le code séquentiel et le code noyau CUDA.
Figure 5 - Comparatif d’exécution entre le code séquentiel sur CPU et le code parallèle sur GPU
À gauche le code séquentiel pris en charge par un thread CPU, à droite le code noyau pris en charge par 1 048 576 threads GPU. Dans le code noyau, les deux boucles externes ont disparu, mais les variables globales blockIdx.y * blockDim.y + threadIdx.y et blockIdx.x * blockDim.x + threadIdx.x, sont renseignées par le moteur d’exécution pour chaque instance de thread CUDA utilisée.
Enfin, nous le signalions dans l’étape 3, nous devons protéger notre code d’itérations excédentaires en filtrant les indices avant de lancer un calcul.
if (j < size && i < size)
{
Ceci peut vous surpendre, mais le coût d’un thread GPU est si faible, qu’un petit excès d’itération ne change pas grand-chose sur le plan des performances.
Étape 5 - Retour du mode GPU/noyau vers le mode CPU/host
L’enchainement des différentes étapes réalisées pour lancer une fonction noyau n’est pas vraiment correct. En effet, les appels noyaux, sont par défaut asynchrones et de ce fait le retour d’un appel noyau est immédiat. Ce mode d’appel est parfaitement justifié, car il permet d’enchainer plusieurs appels noyaux depuis le code CPU/host sans attendre l’exécution effective des calculs. En d’autres mots, le code GPU s’exécute parallèlement sans aucune contrainte de synchronisation.
Étape 6 - Copier les résultats en mode CPU/Host vers la matrice C
Dans l’étape précédente, nous avons précisé que le code noyau s’exécutait de manière asynchrone, mais alors que se passe-t-il, lorsque la fonction cudaMemcpy() est appelée ?
cudaMemcpy(C, deviceC, sizeAlloc, cudaMemcpyDeviceToHost);
Copier la matrice deviceC en provenance du noyau vers la matrice C allouée côté host
La réponse est simple, l’appel à la fonction cudaMemcpy() provoque un point de synchronisation avec le code GPU. Le drapeau cudaMemcpyDeviceToHost assure la copie des octets en provenance de la fonction kernel vers la mémoire host.
Étape 7 - Détruire la mémoire des buffers utilisés par la fonction noyau
La mémoire allouée dans la mémoire globale GPU doit être libérée sous peine de créer une fuite mémoire GPU.
cudaFree(deviceA);
cudaFree(deviceB);
cudaFree(deviceC);
}
Libérer toutes les matrices noyau en fin de traitement
La fonction cudaFree() reprend le contrat de la fonction free() en langage C.
Si vous disposez de tous les codes présentés au sein d’un fichier .cu dans un projet C++ qui respecte le paramétrage présenté dans l’article précédent, alors vous devriez obtenir un binaire opérationnel pour réaliser une exécution.
Exécution du produit matriciel
Après avoir compilé notre code en mode Release en 64bit, nous sommes enfin prêts pour lancer notre produit matriciel avec le Framework CUDA C. Nous pouvons lancer à la fois notre traitement CUDA, que nous appellerons CUDA Simple, précédé des traitements avec un seul cœur et un traitement sur quatre cœurs via la méthode PPL parallel_for présentée dans le premier article.
Figure 6 - Exécution de nos calculs matriciels
Les résultats sont très encourageants, la technologie CUDA C améliore considérablement les performances de notre calcul matriciel.
Figure 7 - Comparatif du calcul matriciel exprimé en millisecondes
Les performances du code CUDA sont extrêmement bonnes vis-à-vis du code parallel_for sur 4 cœurs et incomparablement plus rapide que l’exécution sur un seul cœur. Si vous étiez encore sceptiques sur l’intérêt des GPU, je crois qu'il n’y a plus de doute possible lorsque votre algorithme est massivement orienté données.
Codes sources
L’initiation à la programmation GPU est un sujet délicat qui doit être réalisé avec des codes sources afin de mieux comprendre « comment ça marche ». Je vous invite à télécharger les sources qui m’ont permis d’écrire cet article. Le fichier ZIP contient la solution Visual Studio 2010 « MulMatrix.sln », contenant deux projets : « MulMatrix » et « MulMatrixCuda ».
· Le projet MulMatrix est une simple console, il contient que du code CPU, comme le code séquentiel, parallèle avec PPL et l’appel au code CUDA.
· Le projet MulMatrix est un projet type DLL, il contient à la fois du code host et du code device dans le fichier MulMaltrix.cu.
Point d’entrer de la démonstration.
int _tmain(int argc, _TCHAR* argv[])
{
try
{
print_device();
Demo demo;
demo.run(Demo::SERIAL);
demo.run(Demo::PARALLEL_FOR, false);
demo.run(Demo::PARALLEL_FOR);
demo.run(Demo::CUDA, false);
demo.run(Demo::CUDA);
}
catch (std::exception* e)
{
cerr << e->what();
}
return 0;
}
Pour éviter des mesures de performance altérées par des phénomènes d’initialisation (chargement de DLL, Runtime ..), les traitements sont exécutés deux fois. Une première exécution avec un nombre d’itérations réduites est lancée sans afficher de résultats. Dans le code ci-dessous la méthode run permet d’indiquer si vous souhaitez mesurer ou pas le traitement.
En conclusion
Les architectures CPU et GPU sont fortement différentes et les conséquences pour les développeurs sont importantes. Nous avons appris comment implémenter notre programme matriciel sur une architecture CUDA. Nous avons appris à lancer un traitement noyau en exploitant efficacement la topologie CUDA. Nous avons aussi montré les différentes étapes permettant de lancer, d’exécuter puis de terminer une fonction noyau vis-à-vis de son programme hôte. Et enfin, nous avons comparé les exécutions de la version CUDA de notre algorithme vis-à-vis des versions séquentielles et parallèles sur CPU.
A bientôt,
Bruno