FLASH INFORMATIQUE FI



expert GPU, environnements de programmation CUDA


Cet article s’adresse aux lecteurs évoluant dans le monde du calcul GPU . Avec la sortie des nouvelles cartes NVIDIA Fermi - GTX 470/480, le but ici est de faire le point sur les derniers développements autour de CUDA.



This article is intended for people interested in GPU computing. The new NVDIA Fermi-GTX 470/480 cards having arrived, our purpose here is to examine the latest CUDA developments.


Francis LAPIQUE


Introduction

J’invite le lecteur que le sujet intéresse, mais qui n’a pas de connaissance des architectures des cartes NVIDIA et du modèle de programmation CUDA à lire le livre de NVIDIA et de l’Université de l’Illinois  [1] paru en début ce cette d’année ou à se reporter à une introduction du sujet dans le FI8/09  [2].
Le calcul GPU prend une place de plus en plus importante face au calcul CPU. Au moins deux bonnes raisons pour comprendre ce phénomène relativement récent : les GPU offrent, au niveau applicatif, une bande passante mémoire bien supérieure à celle des CPU et livre, à consommation d’énergie égale, plus de performance calcul.
Quelques annonces prises en vrac :

  • Outils CULA : LAPACK sur les GPU CUDA, développé par EM Photonics
  • MAGMA : LAPACK sur les GPU CUDA et les CPU multi-cœurs, développé par le groupe de recherche de Jack Dongarra
  • GPULib : bibliothèque de fonctions mathématiques pour IDL et MATLAB
  • Bibliothèque de GPU VSIPL (Vector Signal Image Processing Library)
  • OpenCurrent : bibliothèque open-source pour les solveurs d’équations PDE (équations différentielles partielles) accélérés par CUDA sur les grilles régulières
  • Compilateur d’auto-parallélisation de PGI pour Fortran et C vers CUDA C
  • Compilateur d’auto-parallélisation de CAPS HMPP pour C et Fortran vers CUDA C
  • Solveurs PDE Lattice Boltzmann (LBM) sur GPU
  • PyCUDA : wrapper CUDA pour Python
  • Thrust : bibliothèque de modèles C++ pour CUDA
  • R/gputools : ensemble de modules optimisés CUDA (...gpuMatMult gpuMi gpuQr gpuSolve gpuSvd gpuSvmPredict ....)
  • Support d’Eclipse pour CUDA
  • NVIDIA Parallel Nsight.

On ne pourra pas aborder dans le cadre de cet article tous les points de cette liste. Débutons par le dernier, la sortie de NVIDIA Parallel Nsight, ensuite nous parlerons de PyCUDA, de CUDA et LAPACK et de Thrust (C++ et CUDA).

NVIDIA Parallel Nsight

C’est sans doute l’événement le plus important, car c’est celui qui va connaître auprès des communautés de développement autour des GPU l’impact le plus grand.
NVIDIA Parallel Nsight (précédemment connu sous le nom de Nexus) est un environnement de développement pour la programmation d’applications GPU intégré dans Visual Studio 2008 (Windows 7, 32 ou 64-bit).
Cet environnement est ce qui se fait de mieux pour tous celles et ceux qui ont ou auront à programmer des GPU dans un contexte de projet graphique ou de calcul haute performance. Nsight Debugger/Analyzer cible le débogage et l’analyse du code source CUDA C/C++ donc plutôt l’aspect calcul, Nsight Graphics Inspector capture en temps réel des frames Direct3D donc plutôt l’aspect graphique. Nous allons voir quelques exemples dans le domaine du calcul.
Dans la terminologie Nsight il faut faire la distinction entre la machine dite Host qui va héberger votre session Visual Studio et la machine dite Target qui va exécuter le code CUDA. Host et Target peuvent se résumer à une seule machine, mais session et exécution ne peuvent cohabiter sur la même carte GPU. La distinction se fait à l’installation en portant l’extension Nsight sur Host et un programme Monitor sur Target. Host et Target communiquent via le réseau (fig-1).

JPEG - 6.5 ko
fig. 1

Que fait Nsight ? Dans une version standard un code kernel.cu est compilé par le programme nvcc de l’environnement de développement CUDA. Dans le contexte Nsight on fait référence à un autre programme C :\Program Files\NVIDIA Nexus 1.0\CUDA Toolkit\v3.0\Win64\CUDA\bin64\nvcc.exe (cas d’une machine 64 bits) qui va préparer un exécutable pour Visual Studio. Le répertoire v3.0 signale l’utilisation de la dernière version de développement CUDA 3.0.
Avant tout chose, vous devez vous assurer de la bonne configuration de votre projet. Pour la check-list complète, il faut se reporter à la rubrique How To : Work With Samples du guide utilisateur. On peut citer deux ou trois petites choses : vérifier les règles pour l’extension .cu, les propriétés du projet, le chemin d’accès de nvcc. Vous pouvez alors lancer la commande Build, vous avez un compte-rendu dans la fenêtre Output

1>matrixMul - 0 error(s), 1 warning(s)
== Rebuild All: 1 succeeded, 0 failed, 0 skipped ==

Une dernière chose avant de commencer, il faut préciser l’adresse de la machine Target. Sous la rubrique Nexus User Properties on trouve une sous-rubrique Connection Name où vous indiquez l’IP de cette machine. Lancez le programme Monitor sur la machine Target, donnez un point d’arrêt dans votre code, vous pouvez alors passer la commande Start CUDA Debugging sous le menu Nexus. Si la communication entre Target et Host se passe bien le programme s’arrête, dans le cadre de cet exemple, sur l’instruction int tx = threadIdx.x (fig. 2). En pointant tx vous affichez sa valeur. Suivant le modèle de programmation CUDA, un kernel s’exécute sur une grille de blocks eux-mêmes comptant un certain nombre de threads.

JPEG - 4.7 ko
fig. 2

Si vous vous intéressez à un index de blocks et threads particuliers, allez dans le sous-menu Cuda Focus Picker et précisez votre choix (fig. 3).

JPEG - 4.1 ko
fig. 4

Prenons l’exemple d’un produit de matrices (voir la version de l’article en ligne pour le code).
Dans cet exemple on travaille sur une grille de blocks de 8 par 5 et de 16 par 16 de threads par blocks. Mettons-nous dans la peau du thread 0,0 du block 1,0. Pour ce faire on va dans le sous-menu Cuda Focus Picker et on précise notre choix. En cliquant l’icône Step Over pour passer aux instructions suivantes nous arrivons à un passage comme AS(ty, tx) = A[a + wA * ty + tx] où la matrice A est découpée en sous-matrices dans les espaces de shared memory    . Ce que nous voudrions savoir : le thread que nous venons de mettre sous surveillance charge-t-il la bonne sous-matrice  ? Rien de plus simple, ouvrons une fenêtre watch1 (Debug->Windows) et glissons l’expression a + wA * ty + tx dans cette fenêtre pour l’évaluer.

Bien, qu’en est-il des contenus mémoire de la matrice A et sous-matrice As en shared memory. On retourne dans le sous-menu Windows pour ouvrir deux fenêtres Memory. Dans l’une, celle du haut, on affiche A dans l’autre As. N’ayant pas atteint l’instruction __syncthreads() nous n’avons qu’une partie de As (fig. 4). On peut vérifier le chargement complet de As quand on exécute l’instruction de synchronisation.

JPEG - 5.5 ko
fig. 4

Il est possible de demander un point d’arrêt sur une adresse. Dans le sous-menu New Breakpoint de Debug choisir l’option New Data Breakpoint et rentrer une adresse comme &C[0] pour cet exemple. Si on demande la poursuite de l’exécution, on atteint notre nouveau point (fig. 5).

JPEG - 6.8 ko
fig. 5

Une interrogation qui revient souvent quand on développe une application GPU/CUDA les accès mémoire sont-ils coalescents  ? Si ce n’est pas le cas, il faut revoir le code car on se dirige vers des problèmes de performance plus ou moins graves. Pour aborder cette question, il faut passer à un autre aspect du débogage, l’analyse d’activité. Choisir à partir du menu Nexus le sous-menu New Analysis Activity.

JPEG - 5.3 ko
fig. 6

On a coché la case Profile CUDA (fig. 6) et demandé la surveillance des accès non coalescents sur la mémoire globale. Lançons le processus, le résultat est un compte-rendu assez détaillé (fig. 7).

JPEG - 4 ko
fig. 7

En cliquant All sous la rubrique Top Kernels by Time nous avons notre renseignement (fig. 8).

JPEG - 28.5 ko

Signalons que l’on trouve également dans ce compte-rendu une vue de l’API CALL (fig. 9).

JPEG - 5.6 ko
fig. 9

Un dernier aspect, la résolution de problème de violation de mémoire, introduisant dans la source un bug. On active Enable CUDA Memory Checker et on lance l’application. L’application bloque sur le problème et nous indique qu’il a lieu avec le thread (8,15,0) du block (0,4,0) (fig. 10).

PNG - 15 ko
fig. 10

Vous pouvez également aborder les problèmes de divergence de branchement à savoir quand les threads d’un même warp explorent des chemins différents du fait de la présence d’un if-else. En demandant la fenêtre CUDA Device Summary on a une vue sur l’activité des Warps ( groupement de 32 threads). Dans cet exemple nous travaillons sur des blocks de 256 threads on retrouve les 8 warps numérotés de 0 à 7. Ce qui est intéressant c’est le mask de 32 bits (1 thread actif 0 thread inactif ) pour connaître le degré de divergence (fig. 11).

JPEG - 5 ko
fig. 11

Nous n’avons abordé que quelques aspects de cet environnment car le format papier n’est pas le format idéal pour présenter ce genre d’outil, mais insistons sur le fait que c’est un outil majeur pour développer et optimiser des applications massivement parallèles sur GPU.
Passons à un autre point de la liste : la communauté Python et les GPU.

PyCUDA

Pour celles et ceux acquis au langage Python, PyCUDA offre un wrapper CUDA. C’est un environnement idéal pour le prototypage d’applications, car vous allez travailler dans un environnement interprété. Les choses se passent ainsi :
Dans un premier temps, comme pour n’importe quel module, on commence par les import pour l’utilisation de PyCUDA (c’est le moment décisif qui valide ou invalide votre installation....)

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule

On va travailler dans cet exemple sur une matrice 4x4 en simple précision. On importe numpy, qui est le package fondamental pour toute application numérique en Python :

import numpy
a = numpy.random.randn(4,4)
a = a.astype(numpy.float32)

On fait l’allocation mémoire sur le GPU et on transfère la matrice sur le device :

a_gpu = cuda.mem_alloc(a.nbytes)
cuda.memcpy_htod(a_gpu, a)

Le kernel est un code CUDA C. On donne ce source au constructeur pycuda.compiler.SourceModule. Ce kernel, qui porte le nom doublify , double chacune des valeurs de la matrice :

mod = SourceModule("""
__global__ void doublify(float *a)
{
int idx = threadIdx.x + threadIdx.y*4;
a[idx] *= 2;
}
""")

Le code est compilé puis chargé sur le device. La méthode get_function permet d’appeler ce kernel par son nom. On lui donne deux arguments : la matrice et la taille de block 4x4 (par défaut la grille est de 1x1).

func = mod.get_function("doublify")
func(a_gpu, block=(4,4,1))

C’est fini on récupère les données sur le host :

a_doubled = numpy.empty_like(a)
cuda.memcpy_dtoh(a_doubled, a_gpu)

Pour la suite il faut se reporter à la page d’ Andreas Klöckner de Brown University qui conduit ce projet très intéressant. Passons du domaine académique au domaine commercial.

CUDA et LAPACK : CULA

EM Photonics propose des versions optimisées de LAPACK     . La version dite Basic est gratuite mais elle est en simple précision et se limite à 6 fonctions (SGESV SGETRF SGEQRF SGELS SGGLSE SGESVD). La version payante dite Premium est un port complet de LAPACK.
Si vous envisagez une migration, il faut se référer au guide de programmation qui signale les différences d’approche comme le passage par valeur/référence, le rangement orienté ligne/colonne.

// CULA
culaStatus s;
s = culaSgesvd('O', 'N', m, m, a, lda, s, u, ldu, vt, ldvt);

// Traditional
char jobu = 'O';
char jobvt = 'N';
int info;
sgesvd(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &info);

La parallélisation sur plusieurs GPU est, pour le moment, manuelle. Voici un squelette de programation d’une factorisation QR    . Dans le programme principal on va :

  • demander le nombre de devices du système
status = culaGetDeviceCount(&numDevices);
numThreads = numDevices;
  • allouer les structures de données qui vont contenir threads et arguments
threads = (ThreadHandle*)malloc(numThreads*sizeof(ThreadHandle));
args = (ThreadArgs*)malloc(numThreads*sizeof(ThreadArgs));
  • lancer les éxécutions des threads

for(i = 0 ; i < numThreads ; ++i) args[i].id = i ; args[i].gpu = i % numDevices ; args[i].ret = 0 ;

pthread_create(&threads[i], NULL, threadRoutine, &args[i]) ;

  • attendre
                pthread_join(threads[i], NULL);

Pour chacun des threads on va :

  • attacher un GPU
status = culaSelectDevice(args->gpu);
Initialiser CULA
status = culaInitialize();

-*appeler culaSgeqrf

status = culaSgeqrf(M, N, A, M, TAU);
  • quitter CULA

culaShutdown() ; args->ret = EXIT_SUCCESS ;

À titre indicatif voici les résultats de benchmark pour un board Intel® DX48BT2/4GB avec cartes GTX 280 :

Initializing CULA...
Initializing MKL...

Benchmarking the following functions:
SGEQRF SGETRF SGELS SGGLSE SGESV SGESVD
-------------------------------------
-- SGEQRF Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 0.52 2.45 4.6977
5120 0.94 4.58 4.8670
6144 1.55 7.83 5.0615
7168 2.37 12.25 5.1783
8192 2.50 18.32 7.3315
-- SGETRF Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 0.27 1.38 5.1559
5120 0.45 2.56 5.6540
6144 0.71 4.20 5.9301
7168 1.05 6.52 6.2190
8192 1.48 9.65 6.5172
-- SGELS Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 0.74 2.64 3.5537
5120 1.28 4.96 3.8793
6144 2.02 8.31 4.1161
7168 3.00 13.15 4.3880
8192 3.30 19.12 5.7886
-- SGGLSE Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 0.82 6.01 7.3081
5120 1.39 10.24 7.3508
6144 2.20 16.01 7.2857
7168 3.19 23.41 7.3314
8192 3.56 33.00 9.2646
-- SGESV Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 0.33 1.40 4.2376
5120 0.53 2.59 4.9100
6144 0.80 4.23 5.3061
7168 1.15 6.55 5.6845
8192 1.60 9.72 6.0610
-- SGESVD Benchmark --
Size CULA (s) MKL (s) Speedup
------ ---------- ---------- ---------
4096 44.61 208.25 4.6678
5120 49.95 350.64 7.0198
6144 76.12 517.99 6.8048
7168 110.90 820.13 7.3950
8192 154.44 1129.61 7.3144

C++ et CUDA : Thrust

Comme ce titre le laisse entendre Thrust est le pendant de PyCUDA pour les programmeurs de C++. Cette bibliothèque à l’image de la STL   fournit un ensemble de classes de conteneurs telles que host_vector ou device_vector et d’algorithmes comme reduce() ou inclusive_scan(). L’installation est des plus simples puisqu’elle se résume à déposer un répertoire d’include. C’est un environnement conseillé pour prototyper des applications complexes si vous avez une bonne pratique de la STL. Passons rapidement en revue quelques aspects caractéristiques de cette bibliothèque. Les appels aux fonctions cudaMalloc et cudaMemcpy ont été masqués. Les codes gagnent en lisibilité comme le montrent ces quelques lignes :

// allocate host vector with two elements
thrust::host_vector<int> h_vec(2);

// copy host vector to device
thrust::device_vector<int> d_vec = h_vec;

// manipulate device values from the host
d_vec[0] = 13;
d_vec[1] = 27;
std::cout << "sum: " << d_vec[0] + d_vec[1] << std::endl;

host_vector désigne les vecteurs stockés sur CPU, device_vector les vecteurs stockés sur GPU.

On trouve également un ensemble d’itérateurs. Citons le zip_iterator qui permet de construire à peu de frais des structures de données qui ne sont pas naturelles dans CUDA comme des n-uplet "tuple"

// initialize vectors
device_vector<int> A(3);
device_vector<char> B(3);
A[0] = 10; A[1] = 20; A[2] = 30;
B[0] = &#8216;x'; B[1] = &#8216;y'; B[2] = &#8216;z';

// create iterator (type omitted)
first = make_zip_iterator(make_tuple(A.begin(), B.begin()));
last = make_zip_iterator(make_tuple(A.end(), B.end()));

first[0] // returns tuple(10, &#8216;x')

On réduit ces données à la valeur du tuple maximum.

// maximum of [first, last]
maximum< tuple<int,char> > binary_op;
reduce(first, last, first[0], binary_op); // returns tuple(30, &#8216;z')

Calculer la norme d’un vecteur en Thrust c’est réduire les éléments d’un vecteur à une somme dont les éléments ont été élevés au carré ce qui s’écrit de façon relativement formelle comme :

sqrt( transform_reduce(A.begin(), A.end(), unary_op, init, binary_op) )

On a défini une fonction square avec une surcharge l’opérateur () :

struct square
{
__host__ __device__
float operator()(float x) { return x * x; }
};

et les arguments de l’algorithme transform_reduce comme étant :

square unary_op;
plus<float> binary_op;
float init = 0;


/* Matrix multiplication: C = A * B.
* Device code.
*/

#ifndef_MATRIXMUL_KERNEL_H_
#define_MATRIXMUL_KERNEL_H_

#include<stdio.h>
#include "matrixMul.h"

#define CHECK_BANK_CONFLICTS 0
#if CHECK_BANK_CONFLICTS
#define AS(i, j) cutilBankChecker(((float*)&As[0][0]), (BLOCK_SIZE * i + j))
#define BS(i, j) cutilBankChecker(((float*)&Bs[0][0]), (BLOCK_SIZE * i + j))
#else
#define AS(i, j) As[i][j]
#define BS(i, j) Bs[i][j]
#endif

////////////////////////////////////////////////////////////////////////////////
//! Matrix multiplication on the device: C = A * B
//! wA is A's width and wB is B's width
////////////////////////////////////////////////////////////////////////////////
__global__ void
matrixMul( float* C, float* A, float* B, int wA, int wB)
{
   // Block index
   int bx = blockIdx.x;
   int by = blockIdx.y;

   // Thread index
   int tx = threadIdx.x;
   int ty = threadIdx.y;

   // Index of the first sub-matrix of A processed by the block
   int aBegin = wA * BLOCK_SIZE * by;

   // Index of the last sub-matrix of A processed by the block
   int aEnd   = aBegin + wA - 1;

   // Step size used to iterate through the sub-matrices of A
   int aStep  = BLOCK_SIZE;

   // Index of the first sub-matrix of B processed by the block
   int bBegin = BLOCK_SIZE * bx;

   // Step size used to iterate through the sub-matrices of B
   int bStep  = BLOCK_SIZE * wB;

   // Csub is used to store the element of the block sub-matrix
   // that is computed by the thread
   float Csub = 0;

   // Loop over all the sub-matrices of A and B
   // required to compute the block sub-matrix
   for (int a = aBegin, b = bBegin;
            a<= aEnd;
            a += aStep, b += bStep) {

       // Declaration of the shared memory array As used to
       // store the sub-matrix of A
       __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];

       // Declaration of the shared memory array Bs used to
       // store the sub-matrix of B
       __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

       // Load the matrices from device memory
       // to shared memory; each thread loads
       // one element of each matrix
       AS(ty, tx) = A[a + wA * ty + tx];
       BS(ty, tx) = B[b + wB * ty + tx];

       // Synchronize to make sure the matrices are loaded
       __syncthreads();

       // Multiply the two matrices together;
       // each thread computes one element
       // of the block sub-matrix
       for (int k = 0; k<  BLOCK_SIZE; ++k)
       {
           Csub += AS(ty, k) * BS(k, tx);
       }
       // Synchronize to make sure that the preceding
       // computation is done before loading two new
       // sub-matrices of A and B in the next iteration
       __syncthreads();
   }

   // Write the block sub-matrix to device memory;
   // each thread writes one element
   int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
   C[c + wB * ty + tx] = Csub;
}

#endif // #ifndef_MATRIXMUL_KERNEL_H_


Conclusion

Ce panorama est loin d’être complet et nous ferons un second point après la conférence GPU Technology Conference (GTC 2010) qui se tiendra fin septembre. La soumission pour cette conférence vient de s’ouvrir et j’invite tous les membres de la communauté GPU de l’EPFL à participer.

Nous sommes en cours d’évaluation de la consommation électrique par GigaFlops calculé sur des différentes cartes NVIDIA en vue d’évaluer si on peut atteindre le 1 Joule/GFlops sur ces architectures.



Glossaire

coalescence
action de fusionner 2 blocs de mémoire libre adjacents.
GPU
(Graphics Processing Unit) : ce type de processeur graphique, issu du monde des jeux vidéos, est utilisé à présent pour le calcul scientifique, apprécié pour sa structure hautement parallèle.
LAPACK
(Linear Algebra PACKage) : bibliothèque logicielle dédiée à la simulation numérique, très utilisée pour tester les performances (benchmark) d’un processeur pour le calcul numérique.
QR
une des méthodes les plus utilisées en algèbre linéaire pour décomposer une matrice.
shared memory
ou mémoire partagée désigne un large bloc de mémoire vive qui est accédé par différentes unités de calcul au sein d’un ordinateur parallèle.
STL
(Standard Template Library) : cette bibliothèque est un ensemble d’objets et de méthodes standards pour le C++ ; le code C++ normalisé pour certains objets classiques, garantissant un code source lisible, facilement réutilisable

[1] Programming Massively Parallel Processors : A Hands-on Approach David B. Kirk, Wen-mei W. Hwu

[2] Les GPU ne sont pas uniquement faits pour les consoles de jeux du FI 2/2009



Cherchez ...

- dans tous les Flash informatique
(entre 1986 et 2001: seulement sur les titres et auteurs)
- par mot-clé

Avertissement

Cette page est un article d'une publication de l'EPFL.
Le contenu et certains liens ne sont peut-être plus d'actualité.

Responsabilité

Les articles n'engagent que leurs auteurs, sauf ceux qui concernent de façon évidente des prestations officielles (sous la responsabilité du DIT ou d'autres entités). Toute reproduction, même partielle, n'est autorisée qu'avec l'accord de la rédaction et des auteurs.


Archives sur clé USB

Le Flash informatique ne paraîtra plus. Le dernier numéro est daté de décembre 2013.

Taguage des articles

Depuis 2010, pour aider le lecteur, les articles sont taggués:
  •   tout public
    que vous soyiez utilisateur occasionnel du PC familial, ou bien simplement propriétaire d'un iPhone, lisez l'article marqué tout public, vous y apprendrez plein de choses qui vous permettront de mieux appréhender ces technologies qui envahissent votre quotidien
  •   public averti
    l'article parle de concepts techniques, mais à la portée de toute personne intéressée par les dessous des nouvelles technologies
  •   expert
    le sujet abordé n'intéresse que peu de lecteurs, mais ceux-là seront ravis d'approfondir un thème, d'en savoir plus sur un nouveau langage.