CUDA approfondi

Thibaut Cuvelier

Dans l'article précédent, vous avez pu vous familiariser avec CUDA. Nous n'avons que survolé la matière, afin de vous donner goût au calcul par GPU. Cette fois-ci, nous allons un peu plus approfondir certains points, tout en restant dans le runtime, sans nous aventurer dans le driver (pour le moment).

Notez que je partirai du fait que vous avez lu et compris l'article précédent !
Introduction à CUDA

N'hésitez pas à commenter cet article !
Commentez

- Titre : CUDA approfondi

- Auteur : Thibaut Cuvelier

- Parution : 22 avril 2009

- Mise à jour : 9 janvier 2011

- Licence : Copyright ® 2009 Thibaut Cuvelier. Aucune reproduction, même partielle, ne peut être faite de ce site et de l'ensemble de son contenu : textes, documents, images, etc. sans l'autorisation expresse de l'auteur. Sinon vous encourez selon la loi jusqu'à trois ans de prison et jusqu'à 300 000 € de dommages et intérêts.

Sommaire

- I. Introduction

- II. Variables intégrées

- III. Types de vecteurs intégrés

- IV. Textures

- IV-A. Références de texture

- IV-B. Déclaration & attributs à la compilation

- IV-C. Référence des attributs déclarés à l'exécution

- IV-D. Comparaison des mémoires des textures

- IV-D-1. Mémoire linéaire et tableaux

- IV-D-2. Mémoire constante et mémoire globale

- IV-E. Fonctions de fetching

- IV-F. Fonctions de configuration

- IV-G. Remarque

- V. La mémoire

- V-A. La mémoire partagée et les conflits

- V-B. Types d'allocation

- V-C. Mémoire linéaire

- V-D. Tableaux CUDA

- V-E. Copie

- V-F. Résumé des fonctions de copie

- V-F-1. Mémoire linéaire vers mémoire linéaire

- V-F-2. Mémoire linéaire vers tableaux CUDA

- V-F-3. Matrices vers matrices

- V-F-4. Matrices vers tableaux CUDA

- V-F-5. Tableaux CUDA vers mémoire linéaire

- V-F-6. Tableaux CUDA vers matrices

- V-F-7. Tableaux CUDA vers tableaux CUDA

- VI. Les flux

- VI-A. Définition

- VI-B. Création

- VI-C. Événements

- VII. CuBLAS

- VII-A. BLAS

- VII-B. Performances

- VII-C. Les trois niveaux

- VII-D. Initialisation

- VII-E. La mémoire

- VII-F. Construction des nom des de fonctions

- VII-G. Construction des variables

- VII-H. BLAS1

- VII-G-1. Ixamax

- VII-G-2. Ixamin

- VII-G-3. asum

- VII-G-4. axpy

- VII-G-5. copy

- VII-G-6. dot

- VII-G-7. dotc

- VII-G-8. dotu

- VII-G-9. nrm2

- VII-G-10. rot

- VII-G-11. rotg

- VII-G-12. scal

- VII-G-13. swap

- VII-I. BLAS2

- VII-H-1. gbmv

- VII-H-2. gemv

- VII-H-3. sger

- VII-H-4. sbmv

- VII-H-5. spmv

- VII-H-6. spr

- VII-H-7. spr2

- VII-H-8. symv

- VII-H-9. syr

- VII-H-10. syr2

- VII-H-11. tbmv

- VII-H-12. tbsv

- VII-H-13. tpmv

- VII-H-14. tpsv

- VII-H-15. trmv

- VII-H-16. trsv

- VII-J. BLAS3

- VII-I-1. gemm

- VII-I-2. symm

- VII-I-3. syrk

- VII-I-4. syr2k

- VII-I-5. trmm

- VII-I-6. trsm

- VII-K. Gestion des erreurs

- VII-L. Exemple

- VII-M. Voir aussi

- VIII. CuFFT

- VIII-A. La transformée de Fourier

- VIII-B. La transformée rapide de Fourier (FFT)

- VIII-C. Performances

- VIII-D. Plans

- VIII-E. Transformation

- VIII-F. Gestion des erreurs

- VIII-G. Exemples

- VIII-H. Voir aussi

- IX. CUTIL

- IX-A. Macros

- IX-A-1. CUT_DEVICE_INIT

- IX-A-2. CUDA_SAFE_CALL et CUDA_SAFE_CALL_NO_SYNC

- IX-A-3. CU_SAFE_CALL et CU_SAFE_CALL_NO_SYNC

- IX-A-4. CUT_SAFE_CALL

- IX-A-5. CUT_CHECK_ERROR

- IX-A-6. CUT_CHECK_ERROR_GL

- IX-A-7. CUFFT_SAFE_CALL

- IX-A-8. CUT_SAFE_MALLOC

- IX-A-9. CUT_CONDITION

- IX-A-10. CUT_BANK_CHECKER

- IX-B. Fonctions

- IX-B-1. Gestion de la mémoire

- IX-B-1-a. cutFree

- IX-B-2. Gestion des fichiers

- IX-B-2-a. cutFindFilePath()

- IX-B-2-b. cutReadFilef() et dérivés

- IX-B-2-c. cutWriteFilef() et dérivés

- IX-B-2-d. cutLoadPGMub et dérivés, cutLoadPPM(4)ub

- IX-B-2-e. cutSavePGMub et dérivés

- IX-B-3. Gestion de la ligne de commande

- IX-B-3-a. cutCheckCmdLineFlag et dérivés

- IX-B-4. Comparaison

- IX-B-4-a. cutComparef et dérivés

- IX-B-4-b. cutComparefe

- IX-B-5. Gestion du temps

- IX-C. Un exemple complet

- X. CUDPP

- X-A. Qu'est-ce que CUDPP ?

- X-B. Compilation

- X-B-1. Visual Studio

- X-B-2. GCC

- X-C. La découverte par l'exemple

- X-D. Décryptage de l'exemple

- X-D-1. Initialisation

- X-D-2. Les plans

- X-D-3. Le calcul

- X-E. Gestion des erreurs

- X-F. Options de l'algorithme

- X-F-1. Algorithme

- X-F-2. Opérateur

- X-F-3. Type

- X-F-4. Configuration

- X-G. Fonctions de plans

- X-H. Fonctions de calcul

- X-H-1. Scan

- X-H-2. Tri

- XI. Thrust

- XI-A. Qu'est-ce que Thrust ?

- XI-B. Un premier exemple : la version

- XI-C. Les vecteurs

- XI-D. Itérateurs et pointeurs

- XI-E. Algorithmes

- XI-E-1. Transformations

- XI-E-2. Réductions

- XI-E-3. Scans

- XI-E-4. Tri

- XII. Fonctions mathématiques

- XIII. Fonctions atomiques

- XIII-A. Arithmétiques

- XIII-B. Bits

- XIV. Fonction de synchronisation

- XV. Interopérabilité avec DirectX

- XVI. Interopérabilité avec OpenGL

- XVII. Gestion des erreurs

- XVIII. Gestion des périphériques

- XIX. Conclusions

- XIX-A. Téléchargements

- XIX-B. Voir aussi

- XIX-C. Remerciements

I - Introduction

Nous avons étudié quelques principes généraux de l'architecture de CUDA, afin de bien le comprendre. Nous avons ensuite étudié ses différentes briques. Nous nous étions arrêtés sur le runtime, composant que nous allons continuer dans cet article. Ce composant définit quelques variables et types, qui seront bientôt vus. Nous n'avons qu'entr'aperçu la manière d'écrire en mémoire : vous verrez qu'il y a aussi moyen de définir des tableaux en 2D et en 3D, entre autres. Quand vous avez appris à appeler un kernel, je ne vous ai pas dévoilé la consistance du quatrième paramètre de configuration à l'exécution : les flux. Cela sera réglé.

Pour faciliter le développement, NVIDIA met à disposition quelques librairies : CuBLAS et CuFFT. Vous apprendrez à les utiliser, uniquement en C. La communauté GPGPU met à notre disposition une autre librairie, CUDPP, que nous allons aussi aborder. Cette librairie en utilise une autre, CUTIL. Thrust, une dernière librairie, est l'équivalent de la STL pour CUDA.

Mais ces librairies sont de gros mastodontes, parfois, et vous devrez peut-être effectuer de plus petites opérations mathématiques ou atomiques, qui s'adaptent à plus de situations. Cependant, ces calculs ne sont pas synchrones : il faut donc les synchroniser avec votre application.

Et n'oublions pas que les GPU sont avant tout prévus pour la 3D : CUDA peut interopérer avec les API majeures du marché, OpenGL, ainsi que DirectX 9 et 10.

De nos jours, les entreprises ont besoin d'applications fiables, qui ne plantent pas à la moindre erreur : CUDA permet aussi de vérifier les retours des fonctions.

Aussi, les systèmes multi-GPU sont de plus en plus fréquents et CUDA ne permet de n'en utiliser qu'un à la fois, en mode « pilote automatique ». Nous verrons qu'il y a moyen de le choisir. Vous pouvez même utiliser plusieurs GPU en même temps ! Cependant, pas dans le même thread.

II - Variables intégrées

Ces variables sont disponibles dans tous les kernels et indiquent les options de configuration lors du lancement du kernel ainsi que la position du thread dans la grille.

Nom de la variable

Type

Utilité

gridDim

dim3

Dimensions de la grille

blockIdx

uint3

Index du bloc dans la grille

blockDim

dim3

Dimensions du bloc

threadIdx

uint3

Index du thread dans le bloc

warpSize

int

Taille du warp

Cependant, il y a quelques restrictions.

- On ne peut demander l'adresse de ces variables ;

- On ne peut leur affecter une nouvelle valeur.

Voici un exemple d'utilisation.

printf("Dimensions de la grille : ( x = %d ; y = %d ; z= %d )", gridDim.x, gridDim.y, gridDim.z);
printf("Dimensions du bloc : ( x = %d ; y = %d ; z= %d )", blockDim.x, blockDim.y, blockDim.z);
printf("Identifiant de ce bloc : ( x = %d ; y = %d ; z = %d )", blockIdx.x, blockIdx.y, blockIdx.z);
printf("Identifiant de ce thread : ( x = %d ; y = %d ; z = %d )", threadIdx.x, threadIdx.y, threadIdx.z);
printf("Taille du warp : %d", warpSize);

III - Types de vecteurs intégrés

Voici la liste exhaustive de ces types.

- char1

- uchar1

- char2

- uchar2

- char3

- uchar3

- char4

- uchar4

- short1

- ushort1

- short2

- ushort2

- short3

- ushort3

- short4

- ushort4

- int1

- uint1

- int2

- uint2

- int3

- uint3

- int4

- uint4,

- long1

- ulong1

- long2

- ulong2

- long3

- ulong3

- long4

- ulong4

- float1

- float2

- float3

- float4

- double2

Ces types sont dérivés des types du C char, unsigned char, short, unsigned short, int, unsigned int, long, unsigned long, float et double, comme indiqué par leur nom.

Il s'agit de structures. Ils représentent des vecteurs de une à quatre dimensions, accessibles via X, Y, Z et W, dans l'ordre.

On les initialise à l'aide des fonctions TYPE make_TYPE (TYPE x, TYPE y, TYPE z, TYPE w), comme le montre ce court exemple.

float1 floatVectorOne    =    make_float1(1.0);
float2 floatVectorTwo    =    make_float2(1.0, 4.2);
float3 floatVectorThree  =    make_float3(1.0, 4.2, 1475.41742);
float4 floatVectorFour   =    make_float4(1.0, 4.2, 1475.41742, 0.000000000004105);

On peut accéder aux composantes de ces vecteurs très simplement.

floatVectorOne.x;

floatVectorTwo.x;
floatVectorTwo.y;

floatVectorThree.x;
floatVectorThree.y;
floatVectorThree.z;

floatVectorFour.x;
floatVectorFour.y;
floatVectorFour.z;
floatVectorFour.w;

Il existe aussi le type dim3, qui sert à spécifier les dimensions (d'une grille, d'un bloc ...). Il s'agit d'un vecteur à trois dimensions d'entiers. Toutes les composantes non initialisées valent 1.

Nous avons appris dans l'article précédent à les utiliser et à les déclarer. Ce petit bout de code vous rappellera les méthodes à employer.

dim3 dim     (5);
dim3 dimgrid (1, 1);
dim3 dimblock(5, 24, 240);

IV - Textures

Les textures ne sont pas des types semblables aux autres. Elles se déclarent à la manière d'un kernel, mais, avant d'entrer dans ces subtilités, il faut déjà savoir ce qu'est une texture, dans le langage CUDA.

CUDA supporte une partie du matériel de texturing des GPU, utilisé pour les opérations sur la mémoire des textures. Cette mémoire est lue par le kernel quand une fonction de texture fetching (littéralement : aller chercher une texture) est appelée.

Dans cette section, quand je parlerai de tableaux, il s'agira de tableaux CUDA. Vous verrez l'importance de ceci à la prochaine section, concernant la mémoire.

IV-A - Références de texture

Le premier paramètre passé à une telle fonction s'appelle une référence de texture. Une référence de texture précise quels endroits de la mémoire seront utilisés pour cette texture. On doit la lier à une région de la mémoire, la texture, avant de pouvoir l'utiliser. Plusieurs références peuvent pointer sur une même texture, ou sur des textures qui se superposent.

Une référence possède plusieurs attributs. Parmi ceux-ci, la dimensionnalité, qui spécifie la manière d'accéder à la texture : via un tableau à une dimension (avec une seule coordonnée de texture) ; à deux dimensions (avec deux coordonnées de texture) ; ou à trois dimensions (avec trois coordonnées de texture). Les éléments d'une texture sont appelés texels, des éléments de texture (texture elements).

D'autres attributs sont les types des données d'entrée et de sortie, la manière d'interpréter les coordonnées d'entrée, entre autres.

IV-B - Déclaration & attributs à la compilation

Quelques attributs doivent être connus à la compilation et ne peuvent donc pas être changés plus tard. Ces attributs sont définis dès la déclaration de la référence de texture.

texture <Type, Dim, ReadMode> refDeTex;

Type y représente le type de données retourné lors du fetch. Cet attribut est limité aux entiers, aux flottants de simple précision et à tous les types de vecteurs, décrits plus haut.

Dim est la dimensionnalité, le nombre de dimensions, de la texture, qui peut valoir 1, 2 ou 3. Ce paramètre est optionnel et vaut par défaut 1.

ReadMode ne peut prendre que deux valeurs : cudaReadModeNormalizedFloat ou cudaReadModeElementType.

Si Type est un entier de 16 ou de 8 bits et que ReadMode vaut cudaReadModeNormalizedFloat, alors la valeur est, en vérité, retournée en tant que flottant. Toute la plage de valeurs de l'entier est reportée dans l'intervalle [ 0.0 ; 1.0 ] pour un non-signé, dans [ - 1.0 ; + 1.0 ] pour un signé. Par exemple, 0xFF sera lu comme 1.

Si ReadMode vaut cudaReadModeElementType, aucune conversion n'est effectuée.

Il s'agit d'un paramètre optionnel, dont la valeur par défaut est cudaReadModeElementType.

IV-C - Référence des attributs déclarés à l'exécution

Les autres attributs d'une référence sont mutables et peuvent sans problème être changés à l'exécution. Ils spécifient si les coordonnées de la texture sont normalisées ou non, le mode d'adressage et le filtrage.

Par défaut, les textures sont référencées avec des coordonnées flottantes dans l'intervalle [ 0 ; N ], où N est la taille de la texture dans la dimension correspondant à la coordonnée.

Par exemple, une texture de taille 64 x 32 sera référencée avec des coordonnées dans les intervalles [ 0 ; 63 ] et [ 0 ; 31 ] pour les coordonnées x et y.

Des coordonnées normalisées reportent ces intervalles dans l'intervalle [ 0.0 ; 1.0 ]. Cette normalisation convient parfaitement à certaines applications, s'il est nécessaire que les coordonnées soient indépendantes de la taille de la texture, par exemple.

Le mode d'adressage définit ce qui arrive lorsqu'un élément hors dimensions est demandé. Les valeurs en dessous de 0 sont mises à 0 et les variables plus grandes que N sont mises à N-1, dans le cas de coordonnées non normalisées. Dans le cas de coordonnées normalisées, les coordonnées sont ramenées dans l'intervalle [ 0.0 ; 1.0 ]. Dans ce dernier cas, il existe aussi un mode wrap, qui n'utilise que la partie fractionnaire de la coordonnée. 1.25 devient 0.25 ; -1.25 devient 0.75.

Le filtrage de textures linéaires ne peut être effectué que sur des textures qui retournent des flottants. Il effectue des interpolations de basse précision entre les texels proches. Quand ce mode est activé, les texels proches de la cellule recherchée sont lus et la valeur retournée est interpolée, sur base de l'espace entre la valeur recherchée et les données. Une interpolation simple est effectuée pour des textures à une dimension, une interpolation bilinéaire est effectuée pour des textures à deux dimensions.

IV-D - Comparaison des mémoires des textures

IV-D-1 - Mémoire linéaire et tableaux

Une texture peut être stockée en mémoire linéaire ou sous forme de tableaux. La première méthode possède quelques désavantages.

- Dimensionnalité forcée à 1 ;

- Pas de support de filtrage ;

- Adressage uniquement non normalisé ;

- Pas de mode d'adressage : toute valeur hors intervalle est ramenée à 0.

Le matériel met en vigueur une politique d'alignement sur les adresses des textures. Pour rendre ceci plus transparent pour les programmeurs, les fonctions de bind des références renvoient un offset à appliquer aux recherches pour lire la mémoire désirée. Les pointeurs de base retournés par les fonctions d'allocation de CUDA se conforment à cette contrainte, ce qui fait qu'il n'est pas obligatoire de passer l'offset lorsqu'on les utilise.

IV-D-2 - Mémoire constante et mémoire globale

La lecture sur le périphérique avec le principe des textures présente bien des avantages en comparaison des mémoires globale ou constante.

- Elles sont en cache (ce qui améliore fortement les performances si les données sont disponibles) ;

- Elles ne sont pas sujettes aux contraintes sur l'accès à la mémoire que les deux autres doivent respecter pour de bonnes performances ;

- La latence est mieux cachée, améliorant les performances des applications qui lisent aléatoirement la mémoire ;

- Les données d'entrée en entiers sur 8 ou 16 bits peuvent être converties en flottants sur 32 bits de l'intervalle [ 0.0 ; 1.0 ] ou [ -1.0 ; 1.0 ].

Si, en plus, la texture est stockée sous forme de tableau, le matériel fournit d'autres capacités, qui peuvent être utiles pour certaines applications, spécialement dans le domaine du traitement d'images.

Capacité

Utilité

Problème

Filtrage

Interpolation rapide, mais peu précise, de texels

Valide uniquement si la référence de texture renvoie un flottant

Textures à coordonnées normalisées

Codage indépendant de la résolution

Modes d'adressage

Gestion automatique des cas de bordure

Ne peut être utilisé qu'avec des coordonnées normalisées

Cependant, dans le même kernel, le cache des textures n'est pas gardé cohérent en fonction des écritures sur la mémoire globale. Ainsi, toute recherche sur une adresse qui a subi une écriture globale dans le même appel de kernel retourne des données indéfinies.

En d'autres termes, un thread peut lire en toute quiétude un endroit de la mémoire si et seulement si cet endroit a été mis à jour par un autre kernel ou une copie mémoire, mais pas si elle l'a été par le même thread ou un autre thread du même appel de kernel.

Ceci n'a d'utilité que lors de la recherche depuis la mémoire linéaire, comme un kernel ne peut pas écrire dans un tableau.

IV-E - Fonctions de fetching

Dans le cas de la mémoire linéaire, nous utilisons les fonctions de la famille tex1Dfetch().

float tex1Dfetch
(
texture < unsigned char, 1, cudaReadModeNormalizedFloat > refDeTex,
int x
);

x est l'ordonnée du flottant dont on recherche la valeur dans la texture dont une référence est refDeTex.

Le filtrage et les modes d'adressage ne sont pas supportés, contrairement aux couples et aux quadruplets.

float4 tex1Dfetch
(
texture < uchar4, 1, cudaReadModeNormalizedFloat > texRef,
int x
);

Par contre, dans le cas de tableaux, nous utiliserons les familles tex1D(), tex2D() et tex3D().

float tex1D
(
texture < float, 1, cudaReadModeNormalizedFloat > texRef,
float x
);

float tex2D
(
texture < float, 2, cudaReadModeNormalizedFloat > texRef,
float x, float y
);

float tex3D
(
texture < float, 3, cudaReadModeNormalizedFloat > texRef,
float x, float y, float z
);

IV-F - Fonctions de configuration

Le type texture, tel que défini par l'API haut niveau, est un dérivé du type textureReference, lui défini par l'API bas niveau.

struct textureReference
{
int                            normalized     ;
enum    cudaTextureFilterMode  filterMode     ;
enum    cudaTextureAddressMode addressMode[3] ;
struct  cudaChannelFormatDesc  channelDesc    ;
}

normalized non-nul signifie que la texture est normalisée ; toute autre valeur signifie la non-normalisation.

filterMode spécifie le mode de filtrage de la texture. S'il s'agit de cudaFilterModePoint, le texel le plus proche sera retourné ; s'il s'agit de cudaFilterModeLinear, une interpolation sera retournée.

addressMode précise le mode d'adressage : cudaAddressModeWrap impose que les valeurs incorrectes soient remises dans l'ensemble des correctes en n'utilisant que la partie fractionnaire ; cudaAddressModeClamp utilise l'autre mode.
Chacun des éléments du tableau correspond à une dimension à laquelle est appliqué le mode d'adressage.

channelDesc désigne le format de la valeur retournée à la lecture de la texture.

struct cudaChannelFormatDesc
{
int x;
int y;
int z;
int w;
enum cudaChannelFormatKind f;
};

Chacun des entiers donne le nombre de bits de chaque composante retournée, f, leur type.

- cudaChannelFormatKindSigned : entier signé ;

- cudaChannelFormatKindUnsigned : entier non signé ;

- cudaChannelFormatKindFloat : flottant.

Tous ces paramètres peuvent donc être spécifiés à l'exécution depuis l'hôte. Ils ne s'appliquent, comme dit précédemment, qu'aux textures liées à un tableau.

Avant qu'un kernel ne puisse utiliser une référence de texture pour la lire, la référence doit être liée à une texture, avec les fonctions cudaBindTexture() ou cudaBindTextureToArray().

Ce code lie une référence à un espace en mémoire linéaire, pointée par devPtr.

//Avec l'API bas niveau
texture <float, 1, cudaReadModeElementType> texRef;
textureReference * texRefPtr;
cudaGetTextureReference(&texRefPtr, "texRef");
cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
cudaBindTexture(0, texRefPtr, devPtr, &channelDesc, size);

//Avec l'API haut niveau
texture<float, 1, cudaReadModeElementType> texRef;
cudaBindTexture(0, texRef, devPtr, size);

Ce code lie une référence à un tableau, cuArray.

//Avec l'API bas niveau
texture <float, 2, cudaReadModeElementType> texRef;
textureReference * texRefPtr;
cudaGetTextureReference(& texRefPtr, "texRef");
cudaChannelFormatDesc channelDesc;
cudaGetChannelDesc(& channelDesc, cuArray);
cudaBindTextureToArray(texRef, cuArray, & channelDesc);

//Avec l'API haut niveau
texture <float, 2, cudaReadModeElementType> texRef;
cudaBindTextureToArray(texRef, cuArray);

Le format précisé lors de la liaison d'une texture à une référence doit correspondre aux paramètres de la déclaration de référence. Sinon, les résultats sont indéfinis.

La fonction cudaUnbindTexture() sert à délier une référence d'une texture.

IV-G - Remarque

Vous ne pouvez pas déclarer une texture dans un .cu et l'utiliser dans un autre .cu, sous peine de problèmes à l'édition des liens.

Cependant, ceci peut être contourné. Il vous suffit d'utiliser des fonctions comme celles-ci pour les récupérer.

// Declaration d'une texture 2D
texture<float, 2, cudaReadModeElementType> dataTest_tex;

// Accesseur depuis une fonction __host__
texture<float, 2, cudaReadModeElementType>
& getTexture()
{
return dataTest_tex;
}

// Accesseur depuis une fonction __device__ ou __global__
static __inline__ __device__
texture<float, 2, cudaReadModeElementType>
& getDeviceTexture()
{
return dataTest_tex;
}

V - La mémoire

V-A - La mémoire partagée et les conflits

Lorsque nous avons commencé à travailler sur la mémoire partagée, vous n'aviez pas encore beaucoup de connaissances en CUDA, vous n'auriez pas pu comprendre cette section. C'est pourquoi elle est placée ici.

Parce qu'elle est sur la carte elle-même, la mémoire partagée est plus rapide que la mémorie locale et que la mémoire globale. En fait, pour tous les threads d'un warp, y accéder est aussi rapide que pour un registre, tant qu'il n'y a pas de conflit de mémoire entre les threads.

Pour obtenir une bande passante assez haute, la mémoire partagée est divisée en modules de même taille, les banques, auxquelles on peut accéder simultanément. Ainsi, toutes les requêtes en lecture ou écriture faites en n adresses différentes qui tombent dans n banques différentes peuvent être effectuées simultanément, ce qui résulte en une bande passante effective n fois plus élevée que la bande passante d'un seul module.

Cependant, si deux adresses tombent dans la même banque, il y a un conflit de banque et les accès doivent être sérialisés. Le matériel divise une requête en mémoire avec conflits automatiquement en autant de requêtes non conflictuelles que nécessaire. Ceci fait décroître la bande passante effective d'un facteur égal au nombre de requêtes en mémoire à effectuer. Si le nombre de requêtes séparées est de n, la requête en mémoire peut causer jusqu'à n conflits.

Pour des performances maximales, il est donc très important de comprendre comment les adresses mémoire sont associées aux banques, pour pouvoir prévoir les requêtes pour minimiser les conflits.

Dans le cas de la mémoire partagée, les banques sont organisées pour que des mots de 32 bits successifs soient assignés à des banques successives. Chaque banque a donc une bande passante de 32 bits par deux cycles d'horloge.

Pour les périphériques actuels, la taille d'un warp est de 32 threads et il y a 16 banques. Une requête en mémoire partagée pour un warp est divisée en deux requêtes : une pour la première moitié du warp et une autre pour l'autre moitié du warp. Ceci ayant pour conséquence qu'il ne peut y avoir de conflit entre les moitiés d'un warp.

Il est courant que chaque thread accède à un mot de 32 bits d'un tableau indexé, grâce à l'identifiant du thread et à une enjambée de .

__shared__ float shared[32];
float data = shared[BaseIndex + s * tid];

Dans ce cas, les threads et accèdent à la même banque quand est un multiple du nombre de banques , quand est un multiple de , ou lui-même. Ceci se traduit par un conflit uniquement si est impair.

Voici quelques requêtes qui ne produisent pas de conflit.

Contrairement à celles-ci.

D'autres cas méritent d'être mentionnés. Par exemple, il y a conflits si un tableau de char est accédé de cette manière, parce que la taille d'un élément du tableau dépasse ou est en retrait par rapport au mot de 32 bits.

__shared__ char shared[32];
char data = shared[BaseIndex + tid];

Dans cet exemple, shared[1], shared[2], shared[3] et shared[4] sont dans la même banque.

Néanmoins, lorsque l'on y accède de cette manière, il n'y a pas de conflit.

char data = shared[BaseIndex + 4 * tid];

Il y a aussi des conflits pour les tableaux de double, vu que la requête en mémoire est séparée en deux requêtes séparées de 32 bits.

__shared__ double shared[32];
double data = shared[BaseIndex + tid];

Une manière de l'éviter est de séparer le type en deux, de cette manière.

__shared__ int shared_lo[32];
__shared__ int shared_hi[32];

double dataIn;
shared_lo[BaseIndex + tid] = __double2loint(dataIn);
shared_hi[BaseIndex + tid] = __double2hiint(dataIn);

double dataOut = __hiloint2double(shared_hi[BaseIndex + tid],
shared_lo[BaseIndex + tid]);

Cependant, ceci n'améliore pas les performances et devrait même les diminuer sur des architectures futures.

Une assignation de structure est compilée en autant de requêtes en mémoire que nécessaire pour chaque membre de la structure.

__shared__ struct type shared[32];
struct type data = shared[BaseIndex + tid];

Ce code peut avoir beaucoup de résultats différents quant au nombre de requêtes nécessaires.

Définition de type

Conflit à la lecture ?

Pourquoi ?

Non

Chaque membre est lu avec un enjambement d'un mot (32 bits)

Non

Chaque membre est lu avec un enjambement d'un mot (32 bits)

Oui

Chaque membre est lu avec un enjambement de plus d'un mot (5 octets, 40 bits)

Finalement, la mémoire partagée dispose aussi d'un mécanisme de diffusion par lequel un mot de 32 bits peut être lu et diffusé à plusieurs threads simultanément, en n'effectuant qu'une seule opération de lecture. Ceci réduit le nombre de conflits quand plusieurs threads d'un demi-warp lisent une même banque, dans le même mot.

V-B - Types d'allocation

La mémoire, sous CUDA, peut être allouée, soit linéairement, soit sous forme de tableau.

La mémoire linéaire existe dans un espace à adresses de 32 bits. Ainsi, des entités allouées séparément peuvent se référencer l'une l'autre par des pointeurs, par exemple dans un arbre binaire.

Les tableaux sont une mémoire opaque optimisée pour la recherche dans une texture. Ils possèdent d'une à trois dimensions et sont constitués d'éléments à une, deux ou quatre composantes, qui peuvent être des entiers signés ou non sur 8, 16 ou 32 bits, ou bien des flottants sur 32 bits. Ces tableaux sont uniquement lisibles par des kernels avec des recherches dans des textures et ne peuvent être liés qu'à des références avec le même nombre de composantes.

Ces deux types sont accessibles en lecture et écriture par les fonctions de copie mémoire.

V-C - Mémoire linéaire

Elle est allouée avec les fonctions cudaMalloc() et cudaMallocPitch(), effacée avec cudaFree().

Ce code alloue un tableau de 256 éléments flottants en mémoire linéaire.

float* devPtr;
cudaMalloc( (void * *) & devPtr, 256 * sizeof(float) );

Il existe d'autres fonctions pour l'allocation de mémoire : cudaMallocPitch() et cudaMalloc3D(), pour des tableaux en 2D et 3D. Ces fonctions permettent de s'assurer que les copies s'effectueront le plus vite possible, avec les fonctions appropriées : cudaMemcpy2D() et cudaMemcpy3D().

Ce code montre l'allocation et l'utilisation d'un tableau 2D.

float* devPtr;
int pitch;

cudaMallocPitch( (void**)&devPtr, &pitch, width * sizeof(float), height);

mykernel<<<100, 512>>>(devPtr, pitch);

__global__ void mykernel(float* devPtr, int pitch)
{
for (int r = 0 ; r < height ; ++r)
{
float* row = (float *) ( (char *) devPtr + r * pitch);

for (int c = 0 ; c < width ; ++c)
{
float element = row[c];
}
}
}

Et voici pour un tableau 3D.

cudaPitchedPtr devPitchedPtr;
cudaExtent extent = make_cudaExtent(64, 64, 64);

cudaMalloc3D(&devPitchedPtr, extent);

mykernel<<<100, 512>>>(devPitchedPtr, extent);

__global__ void mykernel(cudaPitchedPtr devPitchedPtr, cudaExtent extent)
{
char * devPtr       = devPitchedPtr.ptr;
size_t pitch        = devPitchedPtr.pitch;
size_t slicePitch    = pitch * extent.height;
for (int z = 0; z < extent.depth; ++z)
{
char * slice = devPtr + z * slicePitch;
for (int y = 0 ; y < extent.height ; ++y)
{
float* row = (float*) (slice + y) * pitch;
for (int x = 0 ; x < extent.width ; ++x)
{
float element = row[x];
}
}
}
}

V-D - Tableaux CUDA

Ils sont alloués grâce aux fonctions cudaMallocArray() et cudaMalloc3DArray(), tandis qu'ils sont libérés grâce à cudaFreeArray(). L'allocation nécessite la création d'un descripteur de format, créé à l'aide de la fonction cudaCreateChannelDesc().

cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
cudaArray * cuArray;
cudaMallocArray(& cuArray, & channelDesc, width, height);

La fonction cudaGetSymbolAddress() est utilisée pour récupérer l'adresse pointant sur la mémoire allouée pour une variable qui réside en mémoire globale. La taille peut être demandée grâce à cudaGetSymbolSize()

V-E - Copie

Cet exemple montre la copie d'un tableau linéaire vers un tableau CUDA.

cudaMemcpy2DToArray
(
cuArray,
0,
0,
devPtr,
pitch,
width * sizeof(float),
height,
cudaMemcpyDeviceToDevice
);

Cet exemple montre la copie d'un tableau hôte vers le périphérique.

float data[256];
int size = sizeof(data);
float* devPtr;
cudaMalloc( (void**) & devPtr, size);
cudaMemcpy(devPtr, data, size, cudaMemcpyHostToDevice);

Et celui-ci, vers la mémoire constante.

__constant__ float constData[256];
float data[256];
cudaMemcpyToSymbol(constData, data, sizeof(data));

V-F - Résumé des fonctions de copie

Cette section est complémentaire : chaque type de fonction de copie reprend les paramètres et les retours des fonctions de copie de la mémoire linéaire vers la mémoire linéaire.

V-F-1 - Mémoire linéaire vers mémoire linéaire

cudaError_t cudaMemcpy
(
void * dst,
const void * src,
size_t count,
enum cudaMemcpyKind kind
)

cudaError_t cudaMemcpyAsync
(
void * dst,
const void * src,
size_t count,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
)

- dst : pointeur vers la destination des données ;

- src : source des données ;

- count : nombre de bytes à copier ;

- kind : direction de copie ;

- stream : flux.

kind doit être une de ces valeurs.

- cudaMemcpyHostToHost : hôte à hôte ;

- cudaMemcpyHostToDevice : hôte à périphérique ;

- cudaMemcpyDeviceToHost : périphérique à hôte ;

- cudaMemcpyDeviceToDevice : périphérique à périphérique (ces copies sont automatiquement asynchrones).

cudaMemcpyAsync est la version asynchrone de la copie. Ceci implique un respect de l'hôte. Ainsi, l'appel peut retourner avant que la copie soit terminée. Ceci ne fonctionne qu'en mémoire paginée verrouillée. Si elle n'est pas verrouillée, la fonction retourne une erreur. La copie peut être associée à un flux en en précisant un.

Ces fonctions peuvent retourner ces valeurs.

- cudaSuccess : succès de l'opération ;

- cudaErrorInvalidValue : valeur non valide ;

- cudaErrorInvalidDevicePointer : pointeur sur le périphérique non valable ;

- cudaErrorInvalidMemcpyDirection : mauvaise direction de copie.

V-F-2 - Mémoire linéaire vers tableaux CUDA

cudaError_t cudaMemcpyToArray
(
struct cudaArray * dstArray,
size_t dstX,
size_t dstY,
const void * src,
size_t count,
enum cudaMemcpyKind kind
)

cudaError_t cudaMemcpyToArrayAsync
(
cudaArray * dstArray,
size_t dstX,
size_t dstY,
const void * src,
size_t count,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
)

- dstArray : pointeur sur le tableau CUDA de destination ;

- dstX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- dstY : ordonnée de l'élément à partir duquel la copie doit être effectuée.

V-F-3 - Matrices vers matrices

Les matrices sont aussi appelées tableaux 2D.

cudaError_t cudaMemcpy2D
(
void * dst,
size_t dpitch,
const void * src,
size_t spitch,
size_t width,
size_t height,
enum cudaMemcpyKind kind
)

cudaError_t cudaMemcpy2DAsync
(
void * dst,
size_t dpitch,
const void * src,
size_t spitch,
size_t width,
size_t height,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
)

- dpitch : largeur de la mémoire en bytes de la matrice de la destination ;

- spitch : largeur de la mémoire en bytes de la matrice de la source ;

- width : longueur d'une ligne ;

- height : nombre de lignes.

Cette fonction peut aussi retourner cette erreur.

- cudaErrorInvalidPitchValue : valeurs de pas plus grandes que le maximum autorisé.

V-F-4 - Matrices vers tableaux CUDA

cudaError_t cudaMemcpy2DToArray
(
struct cudaArray * dstArray,
size_t dstX,
size_t dstY,
const void * src,
size_t spitch,
size_t width,
size_t height,
enum cudaMemcpyKind kind
);

cudaError_t cudaMemcpy2DToArrayAsync
(
struct cudaArray * dstArray,
size_t dstX,
size_t dstY,
const void * src,
size_t spitch,
size_t width,
size_t height,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
);

- dstArray : pointeur sur le tableau CUDA de destination ;

- dstX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- dstY : ordonnée de l'élément à partir duquel la copie doit être effectuée ;

- spitch : largeur de la mémoire en bytes de la matrice de la source ;

- width : longueur d'une ligne ;

- height : nombre de lignes.

Cette fonction peut aussi retourner cette erreur.

- cudaErrorInvalidPitchValue : valeurs de pas plus grandes que le maximum autorisé.

V-F-5 - Tableaux CUDA vers mémoire linéaire

cudaError_t cudaMemcpyFromArray
(
void * dst,
const struct cudaArray* srcArray,
size_t srcX,
size_t srcY,
size_t count,
enum cudaMemcpyKind kind
)

cudaError_t cudaMemcpyFromArrayAsync
(
void * dst,
const struct cudaArray* srcArray,
size_t srcX,
size_t srcY,
size_t count,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
)

- srcArray : pointeur sur le tableau CUDA de destination ;

- srcX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- srcY : ordonnée de l'élément à partir duquel la copie doit être effectuée.

V-F-6 - Tableaux CUDA vers matrices

cudaError_t cudaMemcpy2DFromArray
(
void * dst,
size_t dpitch,
const struct cudaArray* srcArray,
size_t srcX,
size_t srcY,
size_t width,
size_t height,
enum cudaMemcpyKind kind
)

cudaError_t cudaMemcpy2DFromArrayAsync
(
void * dst,
size_t dpitch,
const struct cudaArray* srcArray,
size_t srcX,
size_t srcY,
size_t width,
size_t height,
enum cudaMemcpyKind kind,
cudaStream_t stream = 0
)

- dpitch : largeur de la mémoire en bytes de la matrice de la destination ;

- srcArray : pointeur sur le tableau CUDA de destination ;

- srcX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- srcY : ordonnée de l'élément à partir duquel la copie doit être effectuée ;

- width : longueur d'une ligne ;

- height : nombre de lignes.

V-F-7 - Tableaux CUDA vers tableaux CUDA

cudaError_t cudaMemcpyArrayToArray
(
struct cudaArray* dstArray,
size_t dstX,
size_t dstY,
const struct cudaArray * srcArray,
size_t srcX,
size_t srcY,
size_t count,
enum cudaMemcpyKind kind
)

- dstArray : pointeur sur le tableau CUDA de destination ;

- dstX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- dstY : ordonnée de l'élément à partir duquel la copie doit être effectuée ;

- srcArray : pointeur sur le tableau CUDA de destination ;

- srcX : abscisse de l'élément à partir duquel la copie doit être effectuée ;

- srcY : ordonnée de l'élément à partir duquel la copie doit être effectuée ;

Comme vous pouvez le voir, il n'existe pas de version asynchrone de cette fonction !

VI - Les flux

VI-A - Définition

Un flux est une source de données, comme un fichier, la mémoire ou le réseau. Ceci est une définition très générale, qui s'applique à tous les domaines. Dans CUDA, il ne s'agit pas exactement de la même chose.

Pour faciliter les exécutions concurrentes entre hôte et périphérique, certaines fonctions sont asynchrones, dans le runtime (la partie de CUDA que nous apprenons) : l'application récupère le contrôle avant que les calculs soient entièrement effectués. Ainsi, il arrive fréquemment que plusieurs de ces fonctions soient en cours d'exécution en même temps : elles sont concurrentes.

Ces fonctions sont les suivantes.

- cuLaunchGrid() ;

- cuLaunchGridAsync() ;

- Les fonctions de copie suffixées Async ;

- Les copies du périphérique vers le périphérique ;

- Les appels de kernels avec des fonctions __global__.

Les applications s'occupent de la concurrence via des flux. Un flux est une séquence d'instructions, qui doivent s'exécuter dans un certain ordre. D'un autre côté, des flux peuvent arrêter leur exécution pour un autre flux.

VI-B - Création

Un flux est défini en créant un objet flux et en le spécifiant à un appel de kernel, ou à une copie de mémoire. Par exemple, ce code crée deux flux.

cudaStream_t stream[2];
for (int i = 0 ; i < 2 ; ++i)
cudaStreamCreate(& stream[i] );

Chacun de ces flux est défini, par ce code, comme une séquence d'une copie de l'hôte au périphérique, d'un appel de kernel et d'une copie du périphérique vers l'hôte.

for (int i = 0 ; i < 2 ; ++i)
cudaMemcpyAsync(inputDevPtr + i * size, hostPtr + i * size, size, cudaMemcpyHostToDevice, stream[i]);
for (int i = 0 ; i < 2 ; ++i)
mykernel <<< 100, 512, 0, stream[i] >>>    (outputDevPtr + i * size, inputDevPtr + i * size, size);
for (int i = 0 ; i < 2 ; ++i)
cudaMemcpyAsync(hostPtr + i * size, outputDevPtr + i * size, size, cudaMemcpyDeviceToHost, stream[i]);
cudaThreadSynchronize();

Chaque flux copie sa portion du tableau hostPtr vers inputDevPtr dans la mémoire du périphérique, traite ce dernier tableau et copie le résultat outputDevPtr sur l'hôte, dans la partie correspondante de hostPtr.

cudaThreadSynchronize() est appelée à la fin pour s'assurer que tous les flux ont fini de s'exécuter avant d'essayer d'utiliser leur résultat.

cudaStreamSynchronize() peut être appelée pour synchroniser un flux précis, avec tous les autres flux qui continuent leur exécution normale.

cudaStreamDestroy() est utilisée pour la destruction d'un flux.

cudaStreamQuery() permet de vérifier que toutes les opérations du flux déjà envoyées soient bien effectuées.

Tout appel de kernel ou opération de mémoire sans le paramètre de flux ne commence qu'à la fin des autres opérations. Ces opérations sont affectées au flux 0.

VI-C - Événements

Les événements permettent de vérifier l'état d'avancement du flux.

Voici la création de deux événements.

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

Et voici une autre version du code précédent, qui permet de vérifier son état d'avancement, grâce aux événements.

1. cudaEventRecord(start, 0);

3. for (int i = 0; i < 2; ++i)

4. cudaMemcpyAsync(inputDev + i * size, inputHost + i * size, size, cudaMemcpyHostToDevice, stream[i]);

5. for (int i = 0; i < 2; ++i)

6. mykernel <<< 100, 512, 0, stream[i] >>> (outputDev + i * size, inputDev + i * size, size);

7. for (int i = 0; i < 2; ++i)

8. cudaMemcpyAsync(outputHost + i * size, outputDev + i * size, size, cudaMemcpyDeviceToHost, stream[i]);

10. cudaEventRecord(stop, 0);

11. cudaEventSynchronize(stop);

12. float elapsedTime;

13. cudaEventElapsedTime(&elapsedTime, start, stop);

14. cudaEventDestroy(start);

15. cudaEventDestroy(stop);

Ce code est assez explicite : on précise qu'un événement a lieu (ligne 1), puis on exécute les instructions (lignes 3 à 8), on enregistre qu'un autre événement a lieu (ligne 10), on attend qu'il soit émis (ligne 11), on calcule le temps écoulé entre les deux événements (lignes 12 et 13), puis on supprime les événements (lignes 14 et 15).

VII - CuBLAS

VII-A - BLAS

CuBLAS est une implémentation de BLAS pour CUDA, qui s'utilise donc sur GPU. Mais qu'est-ce que BLAS ?

Il s'agit d'un ensemble de fonctions standardisées, initialement prévues pour le Fortran, qui servent au calcul algébrique linéaire basique, comme les multiplications de matrices ou de vecteurs.

BLAS n'est pas une librairie en tant que telle : il ne s'agit que d'un standard, qui doit encore être implémenté. En voici quelques-unes, parmi les principales.

- refblas, l'implémentation de référence ;

- ACML, optimisée par AMD ;

- Accelerate, optimisée par Apple ;

- MKL, optimisée par Intel ;

- CuBLAS, prévue pour CUDA, que nous allons étudier.

Cette librairie, comme les deux autres, n'est prévue que pour être utilisée conjointement à CUDA ! Elle peut être utilisée en solo, c'est pour cela qu'elle a été prévue, mais elle peut aussi très bien utiliser des emplacement de mémoire déjà alloués par CUDA !

VII-B - Performances

L'implémentation qui est proposée par NVIDIA est, évidemment, optimisée pour les GPU NVIDIA. Ce graphique montre les gains que l'on peut espérer en retirer, face à l'ACML, en utilisant Ecolib, une librairie de calcul spécialement prévue pour les GPU, qui utilise CuBLAS.

Comme vous pouvez voir, dans le cas de multiplications de matrices, les gains peuvent être énormes, s'il y a assez de données en entrée. En effet, pour multiplier de petites matrices, les temps d'envoi des données à la carte graphique et d'initialisation du kernel peuvent être démesurés.

Pourquoi avoir choisi cette librairie pour le comparatif ?
Ceux qui l'ont écrite sont des professionnels, qui passent leur temps à optimiser leurs algorithmes, que ce soit pour le GPU ou pour le CPU. On peut donc dire que ce test est fiable, vu que les méthodes de calcul pour CPU et GPU sont optimisées à un très haut niveau, dans les deux cas.

Cependant, lorsque vous ajoutez le fait que les résultats doivent être transportés du GPU vers le CPU, le tableau devient moins parfait.

VII-C - Les trois niveaux

BLAS est constitué de 3 niveaux de fonctions, chacun s'occupant de certains calculs.

Le premier niveau prend en charge les opérations entre vecteurs, comme le produit scalaire, ou la norme, de cette forme.

Le second niveau prend en charge les opérations entre vecteurs et matrices, de cette forme.

Mais aussi les opérations avec une matrice triangulaire de cette forme.

Le troisième niveau prend en charge les opérations entre matrices, de cette forme.

Mais aussi les opérations avec une matrice triangulaire de cette forme.

C'est ce niveau trois qui est utilisé pour la multiplication de matrices générales (alias GEMM).

VII-D - Initialisation

Cette librairie a besoin d'une initialisation en bonne et due forme, mais celle-ci est extrêmement simple.

cublasInit();

Et c'est fait !

Cependant, il peut arriver qu'il y ait des erreurs à l'initialisation. Elles sont renvoyées directement par la fonction. Il n'y a que deux retours possibles à ce stade.

- CUBLAS_STATUS_ALLOC_FAILED : erreur à l'allocation des ressources ;

- CUBLAS_STATUS_SUCCESS : succès.

Et quand la librairie ne sera plus utilisée, il est préférable de l'éteindre.

cublasShutdown();

Ceci peut aussi amener son lot d'erreurs.

- CUBLAS_STATUS_NOT_INITIALIZED : les ressources n'ont pas encore été allouées ;

- CUBLAS_STATUS_SUCCESS : succès.

Néanmoins, la récupération des erreurs des autres fonctions de CUBLAS n'est pas aussi facile : les fonctions ne retournent rien elles-mêmes. Il faut faire appel à cublasGetError() pour récupérer l'erreur. Dès qu'une erreur est lue, l'état est remis à CUBLAS_STATUS_SUCCESS, c'est-à-dire succès.

VII-E - La mémoire

CUBLAS peut s'interfacer avec CUDA même s'il n'est pas fait pour ça à l'origine. Des blocs mémoires alloués avec CUDA peuvent être utilisés avec CuBLAS.

BLAS fonctionne avec deux types de valeurs : des vecteurs et des matrices. CUBLAS permet donc de travailler avec ces types. Avant de leur mettre une valeur, on commence par allouer un espace en mémoire, puis on le spécialise.

Nous allons ici créer une matrice et la remplir.

// Constantes
#define    N    (2);      // Nombre de lignes et de colonnes de la matrice
#define    N2   (N * N);  // Produit du nombre de lignes et de colonnes

// Allocation de la mémoire pour la matrice sur l'hôte
float *      hote;
hote    = (float *) malloc ( N2 * hote[0] );

// Remplissage de la matrice
hote[0] = 14 ;        // 1 1
hote[1] = 54 ;        // 1 2
hote[2] = 29 ;        // 2 1
hote[3] = 36 ;        // 2 2

// Allocation de la mémoire pour la matrice sur le périphérique
float *    periph = 0;
cublasAlloc(N2, sizeof(periph[0]), (void * *) & periph);

// Copie de la matrice de l'hôte vers le périphérique
cublasSetVector(N2, sizeof(hote[0]), hote[0], 1, periph[0], 1);

// Destruction de la matrice sur l'hôte
free(hote);

// Destruction de la matrice sur le périphérique
cublasFree(periph);

CUBLAS ne s'occupe que de matrices dites "colonnes majeures". C'est pourquoi il faut préciser le nombre de lignes. En connaissant le nombre de lignes, on connaît le nombre de colonnes, vu que l'on a l'ensemble des valeurs de la matrice. Ceci nous permet d'utiliser un simple tableau de flottants : chaque valeur du tableau correspond à une valeur de la matrice.

cublasAlloc() crée donc un tableau en mémoire. Les propriétés de ce tableau sont dictées par les paramètres de la fonction.

1. int : nombre de cases du tableau ;

2. int : taille d'un élément du tableau ;

3. void * * : pointeur vers la mémoire à allouer.

cublasSetVector() spécialise ce tableau en vecteur (qui n'est qu'une forme simplifiée de matrice). Voici ses quelques paramètres.

1. int : nombre d'éléments à copier du vecteur hôte ;

2. int : taille d'un élément de chaque côté ;

3. const void * : vecteur hôte ;

4. int : nombre d'espaces entre les valeurs du vecteur hôte ;

5. void * : vecteur périphérique ;

6. int : nombre d'espaces entre les valeurs du vecteur péripérique.

Quand le vecteur est une matrice, ou une partie d'une matrice, l'incrément du vecteur de 1 permet d'accéder à une colonne, de la dimension majeure de la matrice (son nombre de colonnes ou de lignes), à une ligne.

cublasGetVector() fonctionne de manière analogue : le premier vecteur est celui qui va être copié, le second, l'emplacement où le premier va être stocké.

Pour initialiser une vraie matrice, nous pouvons utiliser la fonction cublasSetMatrix(). Voici son prototype.

cublasSetMatrix
(
int rows,
int cols,
int elemSize,
const void *A,
int lda,
void *B,
int ldb
)

- rows : nombre de colonnes ;

- cols : nombre de lignes ;

- elemSize : taille d'un élément de la matrice ;

- A : la matrice source ;

- lda : dimension de direction de A ;

- B : matrice destination ;

- ldb : dimension de direction de B.

Un petit exemple.

// Allocation de la matrice source
float * a;
int rows;
int cols;
int lda;

rows = 4;
cols = 4;
lda  = cols;

a = (float *) malloc( sizeof(float) * rows * cols);

// Remplissage de la matrice source
// Première colonne
a[0] = 1;
a[1] = 2;
a[2] = 3;
a[3] = 4;
// Deuxième colonne
a[rows    ] = 1;
a[rows + 1] = 1;
a[rows + 2] = 1;
a[rows + 3] = 1;
// Troisième colonne
a[rows * 2    ] = 3;
a[rows * 2 + 1] = 4;
a[rows * 2 + 2] = 5;
a[rows * 2 + 3] = 6;
// Quatrième colonne
a[rows * 3    ] = 5;
a[rows * 3 + 1] = 6;
a[rows * 3 + 2] = 7;
a[rows * 3 + 3] = 8;

// Allocation de la matrice destination
float *    b = 0;
cublasAlloc(rows * cols, sizeof(b[0]), (void * *) & b);

// Copie de la matrice
cublasSetMatrix (rows, cols, sizeof(float), (const void *) a, lda, (void *) b, lda) ;

// Suppression des matrices
cublasFree(b);
free(a);

VII-F - Construction des nom des de fonctions

L'implémentation de BLAS a été effectuée à deux niveaux de précision : la précision simple (FP32) et la double précision (FP64).

Il existe aussi des séries de fonctions pour les nombres complexes.

Toutes les combinaisons n'ont pas été implémentées complètement : les fonctions complexes de BLAS1 à double précision et les fonctions complexes de BLAS2 ne le sont pas encore tout à fait. Toutes les autres le sont. Pour plus de précision, cet article précise toutes les combinaisons existantes sous CUDA 2.2.

Ceci a pour conséquence que toutes les fonctions de chaque niveau existent en 4 versions. Voici les préfixes utilisés.

- S : réel, simple précision ;

- C : complexe, simple précision ;

- D : réel, double précision ;

- Z : complexe, double précision.

Ensuite, une série de suffixes, propres à chaque niveau, complètent le nom.

Une fonction de CuBLAS aura donc cet aspect :
cublas + [SDCZ] + suffixe

Par exemple, la fonction du premier niveau qui s'occupe de réels en simple précision et qui effectue un produit scalaire est celle-ci : cublasSdot(). S spécifie le type des vecteurs, dot, l'opération à effectuer.

Il est à noter deux exceptions : les suffixes Ixamax et Ixamin, où le préfixe se met à la place du x.

VII-G - Construction des variables

CuBLAS utilise plusieurs types, deux ne sont pas issus du standard du C : cuComplex et cuDoubleComplex. Vous pouvez les créer grâce à deux fonctions.

cuFloatComplex  make_cuFloatComplex  ( float  x, float  y ) ;
cuDoubleComplex make_cuDoubleComplex ( double x, double y ) ;

Le premier nombre, x, est la partie réelle du complexe ; y, sa partie imaginaire.

Les autres types utilisés sont hérités du standard.

VII-H - BLAS1

VII-G-1 - Ixamax

int cublasIsamax ( int n, const float     * x, int incx ) ;
int cublasIdamax ( int n, const double    * x, int incx ) ;
int cublasIcamax ( int n, const cuComplex * x, int incx ) ;

Trouve l'index de l'élément à la plus grande valeur absolue du vecteur x, qui possède n éléments. Chaque élément est séparé d'un autre de incx.

VII-G-2 - Ixamin

int cublasIsamin ( int n, const float     * x, int incx ) ;
int cublasIdamin ( int n, const double    * x, int incx ) ;
int cublasIcamin ( int n, const cuComplex * x, int incx ) ;

Trouve l'index de l'élément à la plus petite valeur absolue du vecteur x, qui possède n éléments. Chaque élément est séparé d'un autre de incx.

VII-G-3 - asum

float cublasSasum  ( int n, const float    * x, int incx ) ;
float cublasDasum  ( int n, const double   * x, int incx ) ;
float cublasScasum ( int n, const cuDouble * x, int incx ) ;

Retourne la somme des valeurs absolues des éléments du vecteur x, qui possède n éléments. Chaque élément est séparé d'un autre de incx.

VII-G-4 - axpy

void cublasSaxpy ( int n, float     alpha, const float     * x, int incx, float     * y, int incy ) ;
void cublasDaxpy ( int n, double    alpha, const double    * x, int incx, double    * y, int incy ) ;
void cublasCaxpy ( int n, cuComplex alpha, const cuComplex * x, int incx, cuComplex * y, int incy ) ;

Met, dans y, le résultat de l'opération , où est un nombre, et sont des vecteurs.

VII-G-5 - copy

void cublasScopy ( int n, const float     * x, int incx, const float    * y, int incy ) ;
void cublasDcopy ( int n, const double    * x, int incx, const double   * y, int incy ) ;
void cublasCcopy ( int n, const cuComplex * x, int incx, const cuDouble * y, int incy ) ;

Copie le vecteur dans le vecteur .

VII-G-6 - dot

float cublasSdot ( int n, const float  * x, int incx, const float    * y, int incy ) ;
float cublasDdot ( int n, const double * x, int incx, const double   * y, int incy ) ;

Calcule le produit scalaire de et de . En cas d'erreur, cette fonction retourne 0.0f.

VII-G-7 - dotc

cuComplex cublasCdotc ( int n, cuComplex float * x, int incx, const float * y, int incy ) ;

Retourne le produit scalaire des vecteurs et , le premier étant conjugué.

VII-G-8 - dotu

cuComplex cublasCdotu ( int n, cuComplex float * x, int incx, const float * y, int incy ) ;
cuComplex cublasZdotu ( int n, cuComplex float * x, int incx, const float * y, int incy ) ;

Retourne le produit scalaire des vecteurs et .

VII-G-9 - nrm2

float cublasSnrm2  ( int n, const float     * x, int incx ) ;
float cublasDnrm2  ( int n, const double    * x, int incx ) ;
float cublasScnrm2 ( int n, const cuComplex * x, int incx ) ;

Retourne la norme euclidienne du vecteur .

VII-G-10 - rot

void cublasSrot ( int n, float     * x, int incx, float     * y, int incy, float sc, float ss ) ;
void cublasDrot ( int n, double    * x, int incx, double    * y, int incy, float sc, float ss ) ;
void cublasCrot ( int n, cuComplex * x, int incx, cuComplex * y, int incy, float sc, float ss ) ;

Applique une rotation de matrice .

Multiplie la matrice avec la matrice , qui possède n colonnes.

VII-G-11 - rotg

void cublasSrotg ( float     * sa, float     * sb, float * sc, float * ss ) ;
void cublasDrotg ( double    * sa, double    * sb, float * sc, float * ss ) ;
void cublasCrotg ( cuComplex * sa, cuComplex * sb, float * sc, float * ss ) ;

Calcule la matrice de transformation de Givens telle que et que , à partir de sa et sb.

Ceci oblige à calculer , qui sert à calculer les valeurs cherchées de cette manière : .

VII-G-12 - scal

void cublasSscal  ( int n, float           alpha, float           * x, int incx ) ;
void cublasDscal  ( int n, double          alpha, double          * x, int incx ) ;
void cublasCsscal ( int n, cuComplex       alpha, cuComplex       * x, int incx ) ;
void cublasZscal  ( int n, cuDoubleComplex alpha, cuDoubleComplex * x, int incx ) ;

Multiplie le vecteur par la constante .

VII-G-13 - swap

void cublasSrotg ( int n, float     * x, int incx, float * y, int incy ) ;
void cublasDrotg ( int n, double    * x, int incx, float * y, int incy ) ;
void cublasCrotg ( int n, cuComplex * x, int incx, float * y, int incy ) ;

VII-I - BLAS2

VII-H-1 - gbmv

void cublasSgbmv ( char trans, int m, int n, int kl, int ku, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;

Effectue une opération de la forme .

trans doit être une de ces valeurs.

- N ou n : la matrice ne sera pas modifiée avant calcul ;

- T ou t : la transposée de la matrice sera utilisée pour les calculs ;

- C ou c : la transposée de la matrice sera utilisée pour les calculs.

m correspond au nombre de lignes de , n, au nombre de colonnes. Ils doivent être au minimum nuls.

kl correspond au nombre de subdiagonales de la matrice , ku, au nombre de superdiagonales de la matrice .

Diagonale de la matrice (ici, la matrice identité)

Subdiagonale de la matrice

Superdiagonale de la matrice

La matrice a une dimension directrice lda, qui est son nombre de lignes. lda doit au moins être égal à .

est le vecteur qui servira à la multiplication, qui doit avoir au moins dimensions si la matrice n'est pas modifée, ou si sa transposée est utilisée.

est le vecteur qui servira à la multiplication, qui doit avoir au moins dimensions si la matrice n'est pas modifée, ou si sa transposée est utilisée. C'est dans cette variable que sera stocké le résultat de l'opération.

VII-H-2 - gemv

void cublasSgemv ( char trans, int m, int n, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;
void cublasDgemv ( char trans, int m, int n, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;
void cublasZgemv ( char trans, int m, int n, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;

Effectue une opération de la forme . Identique à la fonction précédente, à l'exception des paramètres et .

VII-H-3 - sger

void cublasSger ( int m, int n, float  alpha, const float  * x, int incx, const float  * y,
int incy, float  * A, int lda ) ;
void cublasDger ( int m, int n, double alpha, const double * x, int incx, const double * y,
int incy, double * A, int lda ) ;

Cette fonctione effectue les opérations de cette forme :

VII-H-4 - sbmv

void cublasSsbmv ( char uplo, int n, int k, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;

Effectue les opérations de la forme , où est une matrice de de côté et symmétrique, possédant superdiagonales et subdiagonales. Par exemple, cette matrice peut servir de paramètre, à condition de donner des valeurs à , , et .

uplo doit être une de ces valeurs.

- U ou u : utiliser le triangle supérieur de la matrice (rouge) ;

- L ou l : utiliser le triangle inférieur de la matrice (vert).

VII-H-5 - spmv

void cublasSsbmv ( char uplo, int n, float alpha, const float * AP, const float * x, int incx,
float beta, float * y, int incy ) ;

Effectue les opérations de la forme , identitique à sbmv, à l'exception de la disparition du nombre de superdiagonales.

VII-H-6 - spr

void cublasSspr ( char uplo, int n, float alpha, const float * x, int incx, float * AP ) ;

Effectue les opérations de la forme .

VII-H-7 - spr2

void cublasSspr2 ( char uplo, int n, float alpha, const float * x, int incx,
const float * y, int incy, float * AP ) ;

Effectue les opérations de la forme .

VII-H-8 - symv

void cublasSsymv ( char uplo, int n, float alpha, const float * A, int lda,
const float * x, int incx, float beta, float * y, int incy ) ;

Effectue les opérations de la forme .

VII-H-9 - syr

void cublasSsyr ( char uplo, int n, float alpha, const float * x, int incx,
float * A, int lda ) ;

Effectue les opérations de la forme .

VII-H-10 - syr2

void cublasSsyr2 ( char uplo, int n, float alpha, const float * x, int incx,
const float * x, int incx,  float * A, int lda ) ;

Effectue les opérations de la forme .

VII-H-11 - tbmv

void cublasStbmv ( char uplo, char trans, char diag, int n, const float * A,
int lda, float * x, int incx ) ;

Effectue les opérations de la forme ou bien .

La matrice doit être triangulaire, c'est-à-dire semblable à celle ci-dessous.

Une matrice est dite unitriangulaire si elle est triangulaire et que toutes les valeurs de sa diagonale principale sont 1.

diag permet de spécifier si la matrice est unitriangulaire (U ou u) ou non (N ou n).

VII-H-12 - tbsv

void cublasStbsv ( char uplo, char trans, char diag, int n, int k,
const float * A, int lda, float * x, int incx ) ;

Résout les équations de la forme ou .

VII-H-13 - tpmv

void cublasStpmv ( char uplo, char trans, char diag, int n,
const float * AP, float * x, int incx ) ;

Effectue les opérations de la forme ou bien .

VII-H-14 - tpsv

void cublasStpsv ( char uplo, char trans, char diag, int n,
const float * AP, float * x, int incx ) ;

Résout les équations de la forme ou .

VII-H-15 - trmv

void cublasStrmv ( char uplo, char trans, char diag, int n,
const float * A, int lda, float * x, int incx ) ;

Effectue les opérations de la forme ou bien .

VII-H-16 - trsv

void cublasStrsv ( char uplo, char trans, char diag, int n,
const float * A, int lda, float * x, int incx ) ;

Résout les équations de la forme ou .

VII-J - BLAS3

VII-I-1 - gemm

void cublasSgemm ( char transa, char transb, int m, int n, int k, float alpha,
const float * A, int lda, const float * B, int ldb,
float beta, float * C, int ldc ) ;
void cublasDgemm ( char transa, char transb, int m, int n, int k, double alpha,
const double * A, int lda, const double * B, int ldb,
double beta, double * C, int ldc ) ;
void cublasCgemm ( char transa, char transb, int m, int n, int k, cuComplex alpha,
const cuComplex * A, int lda, const cuComplex * B, int ldb,
cuComplex beta, cuComplex * C, int ldc ) ;
void cublasZgemm ( char transa, char transb, int m, int n, int k, cuDoubleComplex alpha,
const cuDoubleComplex * A, int lda, const cuDoubleComplex * B, int ldb,
cuDoubleComplex beta, cuDoubleComplex * C, int ldc ) ;

Effectue les opérations de la forme .

peut être la matrice , sa transposée , ou bien (dans le cas de matrices complexes) le conjugué de sa transposée .

Les paramètres transa et transb permettent de diriger ce comportement.

- N, n : pas de transformation ;

- T, t : transposée ;

- C, c : conjugué de la transposée.

est une matrice sur .
est une matrice sur .
est une matrice sur .

VII-I-2 - symm

void cublasSsymm ( char side, char uplo, int m, int n, float  alpha, const float  * A, int lda,
const float  * B, int ldb, float  beta, float  * C, int ldc ) ;
void cublasDsymm ( char side, char uplo, int m, int n, double alpha, const double * A, int lda,
const double * B, int ldb, double beta, double * C, int ldc ) ;

Effectue les opérations de la forme (si side vaut L ou l) ou bien (si side vaut R ou r).

et sont des matrices d'ordre fois .

est une matrice symmétrique, de dimension si side vaut L ou l.

VII-I-3 - syrk

void cublasSsyrk ( char uplo, char trans, int n, int k, float alpha,
const float * A, int lda, float beta, float * C, int ldc ) ;
void cublasDsyrk ( char uplo, char trans, int n, int k, double alpha,
const double * A, int lda, double beta, double * C, int ldc ) ;
void cublasZsyrk ( char uplo, char trans, int n, int k, cuDoubleComplex alpha,
const cuDoubleComplex * A, int lda, cuDoubleComplex beta, cuDoubleComplex * C, int ldc ) ;

Effectue les opérations de la forme (si trans vaut N ou n) ou bien (si trans vaut T, t, C, c).

est une matrice carrée d'ordre .

Si trans vaut N ou n, est une matrice d'ordre sur .

Si trans vaut T, t, C ou c, est une matrice d'ordre sur .

VII-I-4 - syr2k

void cublasSsyr2k ( char uplo, char trans, int n, int k, float  alpha,
const float  * A, int lda, const float  * B, int ldb, float  beta, float  * C, int ldc ) ;
void cublasDsyr2k ( char uplo, char trans, int n, int k, double alpha,
const double * A, int lda, const double * B, int ldb, double beta, double * C, int ldc ) ;

Effectue les opérations de la forme (si trans vaut N ou n) ou bien (si trans vaut T, t, C, c).

VII-I-5 - trmm

void cublasStrmm ( char side, char uplo, char transa, char diag, int m, int n, float  alpha,
const float  * A, int lda, const float  * B, int ldb ) ;
void cublasDtrmm ( char side, char uplo, char transa, char diag, int m, int n, double alpha,
const double * A, int lda, const double * B, int ldb ) ;

Effectue les opérations de la forme ou bien .

VII-I-6 - trsm

void cublasStrsm ( char side, char uplo, char transa, char diag, int m, int n,
float  alpha, const float  * A, int lda, float  * B, int ldb ) ;
void cublasDtrsm ( char side, char uplo, char transa, char diag, int m, int n,
double alpha, const double * A, int lda, double * B, int ldb ) ;

Résout les équations de la forme ou bien .

VII-K - Gestion des erreurs

Vous pouvez récupérer la dernière erreur apparue avec la fonction cublasGetError(). Celle-ci va vous renvoyer un des codes suivants (chacun accompagné de sa signification).

Code d'erreur

Signification

CUBLAS_STATUS_SUCCESS

Réussite de l'opération.

CUBLAS_STATUS_NOT_INITIALIZED

CuBLAS n'a pas été initialisé.

CUBLAS_STATUS_ALLOC_FAILED

CuBLAS n'a pas pu allouer de la mémoire (pas assez de mémoire, en général).

CUBLAS_STATUS_INVALID_VALUE

Une variable numérique non supportée a été passée à la fonction.

CUBLAS_STATUS_ARCH_MISMATCH

Le GPU ne supporte pas les opérations en double précision.

CUBLAS_STATUS_MAPPING_ERROR

CuBLAS n'a pas pu accéder au GPU.

CUBLAS_STATUS_EXECUTION_FAILED

CuBLAS n'a pas pu exécuter votre demande correctement sur le GPU.

CUBLAS_STATUS_INTERNAL_ERROR

Erreur interne.

VII-L - Exemple

De l'utilisation de CuBLAS pour la multiplication de deux matrices aléatoires.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "cublas.h"

// Taille des matrices
#define N (275)
#define N2 (N * N);

int main(int argc, char** argv)
{
// Déclaration des variables
cublasStatus status;
float *    h_A;
float *    h_B;
float *    h_C;
float *    d_A        = 0;
float *    d_B        = 0;
float *    d_C        = 0;
float    alpha    = 1.0f;
float    beta    = 0.0f;
int        i;

// Initialisation de CuBLAS
// avec gestion des erreurs
status = cublasInit();
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur lors de l'initialisation de CuBLAS !!!!!\n");
return EXIT_FAILURE;
}

// Allocation de la mémoire pour les matrices sur l'hôte
// avec gestion des erreurs
//h_A
h_A = (float *) malloc ( N2 * sizeof(h_A[0]) );
if (h_A == 0)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice A sur l'hôte !!!!!\n");
return EXIT_FAILURE;
}
//h_B
h_B = (float * ) malloc ( N2 * sizeof(h_B[0]) );
if (h_B == 0)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice B sur l'hôte !!!!!\n");
return EXIT_FAILURE;
}
//h_C
h_C = (float * ) malloc ( N2 * sizeof(h_C[0]) );
if (h_C == 0)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice C sur l'hôte !!!!!\n");
return EXIT_FAILURE;
}

//Remplissage des matrices avec des données aléatoires
for (i = 0 ; i < N2 ; i++)
{
h_A[i] = rand() / (float) RAND_MAX;
h_B[i] = rand() / (float) RAND_MAX;
h_C[i] = rand() / (float) RAND_MAX;
}

// Allocation de la mémoire pour les matrices sur le périphérique
// avec gestion des erreurs
//d_A
status = cublasAlloc (N2, sizeof (d_A[0]), (void **) & d_A);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice A sur le périphérique !!!!!\n");
return EXIT_FAILURE;
}
//d_B
status = cublasAlloc (N2, sizeof (d_B[0]), (void **) & d_B);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice B sur le périphérique !!!!!\n");
return EXIT_FAILURE;
}
//d_C
status = cublasAlloc (N2, sizeof (d_C[0]), (void **) & d_C);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur lors de l'allocation mémoire de la matrice C sur le périphérique !!!!!\n");
return EXIT_FAILURE;
}

// Remplissage des matrices du périphérique avec les données des matrices de l'hôte
// avec gestion des erreurs
status = cublasSetVector(N2, sizeof(h_A[0]), h_A, 1, d_A, 1);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur à l'accès du périphérique lors de l'écriture de la matrice A !!!!!\n");
return EXIT_FAILURE;
}

status = cublasSetVector(N2, sizeof(h_B[0]), h_B, 1, d_B, 1);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur à l'accès du périphérique lors de l'écriture de la matrice A !!!!!\n");
return EXIT_FAILURE;
}

status = cublasSetVector(N2, sizeof(h_C[0]), h_C, 1, d_C, 1);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!!! Erreur à l'accès du périphérique lors de l'écriture de la matrice A !!!!!\n");
return EXIT_FAILURE;
}

// Nettoyage de la mémoire des erreurs
cublasGetError();

// Multiplication des matrices avec CuBLAS
cublasSgemm('n', 'n', N, N, N, alpha, d_A, N, d_B, N, beta, d_C, N);

// Gestion des erreurs
status = cublasGetError();
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!! Erreur lors de l'exécution du kernel CuBLAS !!!!!\n");
return EXIT_FAILURE;
}

// Allocation de la mémoire sur l'hôte pour la récupération du résultat depuis le périphérique
// avec gestion des erreurs
h_C = (float*)malloc(N2 * sizeof(h_C[0]));
if (h_C == 0) {
fprintf (stderr, "!!!!! Erreur lors de l'allocation de la mémoire sur l'hôte pour la matrice C !!!!!\n");
return EXIT_FAILURE;
}

// Lecture du résultat
// avec gestion des erreurs
status = cublasGetVector(N2, sizeof(h_C[0]), d_C, 1, h_C, 1);
if (status != CUBLAS_STATUS_SUCCESS) {
fprintf (stderr, "!!!!! Erreur d'accès au périphérique (lecture de la matrice C) !!!!!\n");
return EXIT_FAILURE;
}

// Libération de la mémoire sur l'hôte
free(h_A);
free(h_B);
free(h_C);

// Libération de la mémoire sur le périphérique
// avec gestion des erreurs
status = cublasFree(d_A);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!! Erreur lors de la libération de l'espace périphérique pour la matrice A !!!!!\n");
return EXIT_FAILURE;
}
status = cublasFree(d_B);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!! Erreur lors de la libération de l'espace périphérique pour la matrice B !!!!!\n");
return EXIT_FAILURE;
}
status = cublasFree(d_C);
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!! Erreur lors de la libération de l'espace périphérique pour la matrice C !!!!!\n");
return EXIT_FAILURE;
}

// Fermeture propre de CuBLAS
status = cublasShutdown();
if (status != CUBLAS_STATUS_SUCCESS)
{
fprintf (stderr, "!!!! Erreur à la fermeture de CuBLAS !!!!!\n");
return EXIT_FAILURE;
}

printf("\n Appuyez sur ENTER pour quitter.\n");
getchar();

return EXIT_SUCCESS;
}

Vous avez pu remarquer que l'entièreté de CuBLAS est contenue dans le fichier cublas.h.

VII-M - Voir aussi

- refblas

- AMD ACML

- Apple Accelerate

- Intel MKL

- Ecolib

- Algèbre linéaire

- BLAS

- Produit scalaire

- Norme

- Norme euclidienne

- Nombre complexe

- Rotation plane

- Rotation de Givens

- Matrices et quaternions

- Matrice

- Ordre d'une matrice

- Matrice identité

- Diagonale d'une matrice

- Multiplication de matrices générales (GEMM)

- Transposée

- Matrice de rotation

- Conjugué

VIII - CuFFT

VIII-A - La transformée de Fourier

Il s'agit d'un algorithme de décomposition d'une séquence de valeurs en composants de différentes fréquences.

En pratique, cela peut servir à décomposer un signal audio pour obtenir la répartition des fréquences en fonction du temps, comme sur votre chaîne Hi-Fi, qui peut décomposer le signal transmis par les baffles en plages de fréquences : 40 Hz, 200 Hz, 1 kHz, 5kHz et 12 kHz, en général. Mais aussi pour votre électrocardiogramme !

La transformée permet d'obtenir la répartition des fréquences qui caractérisent toutes les ondes qui composent l'onde totale. Le résultat s'appelle un diffractogramme.

Cet exemple se limite aux éléments à une seule dimension. Or, notre univers est constitué de trois dimensions visibles. Cette transformée peut aussi s'appliquer à des éléments à deux ou trois dimensions.

Un exemple d'élément à deux dimensions ? Une image, tout simplement. La transformée permet de passer de l'espace réel (décrit avec des distances) à l'espace réciproque (décrit avec des fréquences). Les plus fins détails de l'image seront transcrits par les plus hautes fréquences de la réciproque.

Concernant la troisième dimension, cette transformée peut servir à la résolution d'équations différentielles partielles en trois dimensions.

VIII-B - La transformée rapide de Fourier (FFT)

Il s'agit d'un algorithme de calcul de la transformée discrète de Fourier, nettement plus rapide. Cet algorithme est le plus utilisé en traitement numérique de signal.

Cette transformée est un cas particulier des transformées de Fourier, qui translatent une fonction complexe d'une variable réelle en une autre fonction, la représentation du domaine fréquentiel (ou de fréquence), qui décrit quelles sont les fréquences dans la fonction d'origine.

Ces algorithmes (car il y en a plusieurs) sont plus rapides, mais aussi parfois plus précis dans certains cas, d'autres algorithmes permettent une approximation très approchée des valeurs réelles.

CuFFT permet, à l'heure actuelle, des transformations en une, deux ou trois dimensions, de valeurs réelles (flottantes) ou bien complexes. Des transformations à une dimension peuvent être effectuées en parallèle.

Les données à traiter sont limitées quant à leur taille : pour une seule dimension, seuls huit millions d'éléments sont supportés. Pour deux et trois dimensions, les tailles sont limitées à l'intervalle [ 2 ; 16384 ].

VIII-C - Performances

Ici, CuFFT sera comparé à FFTW, une des librairies les plus rapides pour CPU, selon les quelques comparatifs produits avec BenchFFT.

Comparatif entre FFTW et CuFFT

Nous pouvons voir très aisément que le CPU reste en très bonne place, pour de petites quantités de données. Cependant, lorsque l'on commence à avoir des paquets plus importants, le GPU s'impose très vite.

À propos de ce GPU, le graphique montre aussi à quel point les transferts entre CPU et GPU peuvent être lents : avec ces transferts, beaucoup moins de données sont traitées, mais toujours plus qu'un CPU, à partir d'un certain point.

Comparatif entre CuFFT et MKL

Ici, CuFFT est comparé à Intel MKL, une librairie prévue pour les mathématiques (sur processeurs Intel, pour de meilleurs résultats). Ici, côté CPU, on pourrait dire que l'algorithme est très fort optimisé, mais ce n'est pas suffisant pour détrôner le GPU, qui reste le plus rapide.

VIII-D - Plans

CuFFT est basé sur FFTW. Ainsi, les fonctions ne traitent que des pointeurs sur des plans (* cufftHandle).

Mais qu'est-ce qu'un plan ? Un plan est une configuration pour effectuer un calcul, connue à l'avance, ce qui permet d'optimiser un maximum les calculs.

Il existe trois fonctions pour créer des plans : cufftPlan1d(), cufftPlan2d() et cufftPlan3d(), dont voici les prototypes.

cufftResult cufftPlan1d ( cufftHandle * plan, int nx,                 cufftType type, int batch ) ;
cufftResult cufftPlan3d ( cufftHandle * plan, int nx, int ny,         cufftType type            ) ;
cufftResult cufftPlan2d ( cufftHandle * plan, int nx, int ny, int nz, cufftType type            ) ;

Le plan qui sera créé aura un pointeur plan. Il contiendra nx éléments sur l'axe des abscisses, ny sur l'axe des ordonnées et nz sur l'axe des cotes (suivant le nombre de dimensions du plan). Le plan sera transformé batch fois (seulement dans le cas de plans à une dimension). Le type de transformation sera type, qui peut avoir ces quelques valeurs.

- CUFFT_C2C : complexe vers complexe ;

- CUFFT_C2R : complexe vers réel ;

- CUFFT_R2C : réel vers complexe.

Un plan peut être détruit grâce à cufftResult cufftDestroy ( cufftHandle plan ) ;.

Les valeurs des plans peuvent être, soit réelles, soit complexes. Les types associés sont cufftReal et cufftComplex.

VIII-E - Transformation

Il existe trois fonctions permettant d'effectuer la transformation, correspondant aux trois valeurs possibles de type pour les fonctions de création de plan.

cufftResult cufftExecC2C ( cufftHandle plan, cufftComplex * idata, cufftComplex * odata, int direction ) ;
cufftResult cufftExecR2C ( cufftHandle plan, cufftReal    * idata, cufftComplex * odata                ) ;
cufftResult cufftExecC2R ( cufftHandle plan, cufftComplex * idata, cufftReal    * odata                ) ;

idata est un pointeur vers les données à traiter, odata est un pointeur vers l'emplacement où les valeurs traitées seront stockées.

direction spécifie la direction de l'algorithme : soit vers l'avant (CUFFT_FORWARD), soit vers l'arrière (CUFFT_BACKWARD).

VIII-F - Gestion des erreurs

Contrairement à CUDA et à CuBLAS, CuFFT n'utilise pas une fonction spécifique pour récupérer les erreurs. Chaque fonction renvoie un cufftResult, qui peut avoir une de ces valeurs.

- CUFFT_SUCCESS : opération effectuée correctement ;

- CUFFT_INVALID_PLAN : plan passé en argument invalide ;

- CUFFT_ALLOC_FAILED : erreur lors de l'allocation de la mémoire ;

- CUFFT_INVALID_TYPE : type demandé non supporté ;

- CUFFT_INVALID_VALUE : mauvais pointeur mémoire ;

- CUFFT_INTERNAL_ERROR : erreur interne ;

- CUFFT_EXEC_FAILED : erreur d'exécution de la transformée sur le GPU ;

- CUFFT_SETUP_FAILED : erreur d'initialisation de la librairie ;

- CUFFT_SHUTDOWN_FAILED : erreur de fermeture de la librairie ;

- CUFFT_INVALID_SIZE : taille de plan non supportée.

VIII-G - Exemples

Pour une transformation de complexe vers complexe à une dimension, sans gestion des erreurs.

#define NX (256)
#define BATCH (10)
#define NB (NX * BATCH)

int main()
{
cufftHandle plan;
cufftComplex *data;
cudaMalloc( (void **) & data, sizeof (cufftComplex) * NB);

cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH);

cufftExecC2C(plan, data, data, CUFFT_FORWARD);

cufftDestroy(plan);
cudaFree(data);

return EXIT_SUCCESS;
}

Pour une transformation de complexe vers complexe à deux dimensions.

#define NX  (256)
#define NY  (128)
#define NXY (NX * NY)

int main()
{
cufftHandle    plan;
cufftComplex * idata;
cufftComplex * odata;

cudaMalloc( (void **) & idata, sizeof (cufftComplex) * NXY);
cudaMalloc( (void **) & odata, sizeof (cufftComplex) * NXY);

status = cufftPlan2d(& plan, NX,NY, CUFFT_C2C);
if (status != CUFFT_SUCCESS)
return EXIT_FAILURE;

status = cufftExecC2C(plan, idata, odata, CUFFT_FORWARD);
if (status != CUFFT_SUCCESS)
return EXIT_FAILURE;

status cufftDestroy(plan);
if (status != CUFFT_SUCCESS)
return EXIT_FAILURE;

cudaFree(idata);
cudaFree(odata);

return EXIT_SUCCESS;
}

Ces exemples ne disent pas une chose : pour pouvoir utiliser CuFFT, il faut inclure le fichier cufft.h.

VIII-H - Voir aussi

- La transformée de Fourier en traitement d'images

- La transformée de Fourier et la transformée rapide de Fourier

- FFTW

- Utilisation de FFTW

- BenchFFT

IX - CUTIL

CUDA Utility Library.

Afin de bien comprendre CUDPP, il faut connaître, au moins un peu, CUTIL, qui simplifie le développement. Cette librairie est, à l'origine, prévue pour raccourcir le temps de développement des exemples du SDK. Ainsi, elle supporte diverses choses :

- Parser la ligne de commande ;

- Lire et écrire des fichiers binaires et des images aux formats PPM et PGM ;

- Comparer des tableaux de données (pour la comparaison des résultats entre CPU et GPU) ;

- Mesurer le temps ;

- Chercher des erreurs dans le code, par des macros ;

- Rechercher les conflits dans les banques de mémoire partagée.

Comme elle est utilisée dans CUDPP et que sa compilation est fort simple, elle est reléguée à la prochaine partie.

Toutes ces macros et fonctions sont disponibles dans les fichiers cutil.h et cutil_gl_error.h

IX-A - Macros

Ces macros ne remplissent leur office que lorsque l'application est compilée en mode de débogage ! Si elle est compilée en Release, la macro n'aura aucun impact sur les performances.
En cas d'erreur, elles ferment automatiquement l'application, en même temps qu'elles écrivent dans le flux stderr.

IX-A-1 - CUT_DEVICE_INIT

Cette macro ne prend pas de paramètres. Elle initialise le premier périphérique disponible pour CUDA. En cas de compilation en mode d'émulation, cette fonction ne fait rien.

CUT_DEVICE_INIT();

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

There is no device.

Ou bien, selon le cas.

There is no device supporting CUDA.

Cette macro est la seule à ne pas changer en fonction du mode de compilation !

IX-A-2 - CUDA_SAFE_CALL et CUDA_SAFE_CALL_NO_SYNC

Ces macros prennent, comme paramètre, un appel à une fonction du runtime de CUDA. En cas d'erreur, elles écrivent l'erreur dans stderr et ferment proprement l'application.

La première commence par synchroniser les threads, la seconde passe cette étape.

CU_SAFE_CALL( cuFree );
CU_SAFE_CALL_NO_SYNC( cuFree );

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

Cuda runtime error %x in file __FILE__ in line __LINE__.
Cuda runtime error %x in file __FILE__ in line __LINE__.

IX-A-3 - CU_SAFE_CALL et CU_SAFE_CALL_NO_SYNC

Ces macros prennent, comme paramètre, un appel à une fonction du driver de CUDA. En cas d'erreur, elles écrivent l'erreur dans stderr et ferment proprement l'application.

CU_SAFE_CALL( cuFree );
CU_SAFE_CALL_NO_SYNC( cuFree );

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

Cuda driver error %x in file __FILE__ in line __LINE__.
Cuda driver error %x in file __FILE__ in line __LINE__.

IX-A-4 - CUT_SAFE_CALL

Cette macro prend, comme paramètre, un appel à une fonction du driver de CUDA. En cas d'erreur, elles écrivent l'erreur dans stderr et ferment proprement l'application.

CUT_SAFE_CALL( cutFree );

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

Cut error in file __FILE_ in line __LINE__.

IX-A-5 - CUT_CHECK_ERROR

Cette macro prend en paramètre le message d'erreur à afficher dans stderr.

CUT_CHECK_ERROR("'not enough memory'");

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

Cuda error: 'not enough memory' in file __FILE__ in line __LINE__ : %x.

IX-A-6 - CUT_CHECK_ERROR_GL

Cette macro ne prend pas de paramètre et affiche l'erreur OpenGL.

CUT_CHECK_ERROR_GL();

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

GL Error in file __FILE__ in line __LINE__ :
%s

IX-A-7 - CUFFT_SAFE_CALL

Cette macro prend, en paramètre, un appel à une fonction de CuFFT et affiche une erreur, s'il y a lieu.

CUFFT_SAFE_CALL( cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH) );

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

CUFFT error in file __FILE__ in line __LINE__.

IX-A-8 - CUT_SAFE_MALLOC

Cette macro prend en paramètre un appel à malloc() et retourne, le cas échéant, une erreur dans le flux stderr.

CUT_SAFE_MALLOC( cudaMalloc( (void **) & buffer, size)) ;

Résultera, en cas d'erreur, en ceci, dans le flux stderr.

Host malloc failure in file __FILE__ in line __LINE__.

IX-A-9 - CUT_CONDITION

Prend, en paramètre, une condition et ferme l'application si elle n'est pas remplie.

CUT_CONDITION(0 == 0);

IX-A-10 - CUT_BANK_CHECKER

IX-B - Fonctions

IX-B-1 - Gestion de la mémoire

IX-B-1-a - cutFree

Cette fonction détruit le pointeur sur une mémoire allouée avec CUTIL qui lui est passé comme paramètre. Si un pointeur a la valeur NULL, cette fonction doit d'abord être utilisée sur lui avant qu'il ne puisse être utilisé avec une autre fonction.

void cutFree ( void * ptr ) ;

IX-B-2 - Gestion des fichiers

IX-B-2-a - cutFindFilePath()

Trouve le chemin du fichier filename, le chemin vers l'exécutable de votre application étant executablePath. La fonction retourne, soit le chemin vers le fichier, soit 0.

char * cutFindFilePath ( const char * filename, const char * executablePath ) ;

IX-B-2-b - cutReadFilef() et dérivés

Retourne, dans data, pointeur non initialisé, un pointeur vers les données de filename qui ont été lues. len est la longueur de ces données. La fonction retourne CUTTrue si tout s'est bien passé, CUTFalse sinon.

CUTBoolean cutReadFilef  ( const char * filename, float         ** data, unsigned int * len,
bool verbose = false ) ;
CUTBoolean cutReadFiled  ( const char * filename, double        ** data, unsigned int * len,
bool verbose = false ) ;
CUTBoolean cutReadFilei  ( const char * filename, int           ** data, unsigned int * len,
bool verbose = false ) ;
CUTBoolean cutReadFileui ( const char * filename, unsigned int  ** data, unsigned int * len,
bool verbose = false ) ;
CUTBoolean cutReadFileb  ( const char * filename, char          ** data, unsigned int * len,
bool verbose = false ) ;
CUTBoolean cutReadFileub ( const char * filename, unsigned char ** data, unsigned int * len,
bool verbose = false ) ;

IX-B-2-c - cutWriteFilef() et dérivés

Retourne, dans data, pointeur non initialisé, un pointeur vers les données de filename qui ont été lues. len est la longueur de ces données. epsilon sert à préciser la précision. La fonction retourne CUTTrue si tout s'est bien passé, CUTFalse sinon.

CUTBoolean cutWriteFilef  ( const char * filename, const float         * data,
unsigned int len, float         epsilon, bool verbose = false ) ;
CUTBoolean cutWriteFiled  ( const char * filename, const double        * data,
unsigned int len, double        epsilon, bool verbose = false ) ;
CUTBoolean cutWriteFilei  ( const char * filename, const int           * data,
unsigned int len, int           epsilon, bool verbose = false ) ;
CUTBoolean cutWriteFileui ( const char * filename, const unsigned int  * data,
unsigned int len, unsigned int  epsilon, bool verbose = false ) ;
CUTBoolean cutWriteFileb  ( const char * filename, const char          * data,
unsigned int len, char          epsilon, bool verbose = false ) ;
CUTBoolean cutWriteFileub ( const char * filename, const unsigned char * data,
unsigned int len, unsigned char epsilon, bool verbose = false ) ;

IX-B-2-d - cutLoadPGMub et dérivés, cutLoadPPM(4)ub

Lit l'image file dans le pointeur data. Elle possède h pixels de large et w pixels de large.

CUTBoolean cutLoadPGMub  ( const char * file, unsigned char  ** data, unsigned int * w,
unsigned int * h ) ;
CUTBoolean cutLoadPPMub  ( const char * file, unsigned char  ** data, unsigned int * w,
unsigned int * h ) ;
CUTBoolean cutLoadPPM4ub ( const char * file, unsigned char  ** data, unsigned int * w,
unsigned int * h ) ;

CUTBoolean cutLoadPGMi   ( const char * file, unsigned int   ** data, unsigned int * w,
unsigned int * h ) ;

CUTBoolean cutLoadPGMs   ( const char * file, unsigned short ** data, unsigned int * w,
unsigned int * h ) ;

CUTBoolean cutLoadPGMf   ( const char * file, unsigned float ** data, unsigned int * w,
unsigned int * h ) ;

IX-B-2-e - cutSavePGMub et dérivés

Écrit l'image file depuis le pointeur data. Elle possède w pixels de large et h pixels de large.

CUTBoolean cutSavePGMub  ( const char * file, unsigned char  * data, unsigned int w, unsigned int h ) ;
CUTBoolean cutSavePPMub  ( const char * file, unsigned char  * data, unsigned int w, unsigned int h ) ;
CUTBoolean cutSavePPM4ub ( const char * file, unsigned char  * data, unsigned int w, unsigned int h ) ;

CUTBoolean cutSavePGMi   ( const char * file, unsigned int   * data, unsigned int w, unsigned int h ) ;

CUTBoolean cutSavePGMs   ( const char * file, unsigned short * data, unsigned int w, unsigned int h ) ;

CUTBoolean cutSavePGMf   ( const char * file, unsigned float * data, unsigned int w, unsigned int h ) ;

IX-B-3 - Gestion de la ligne de commande

IX-B-3-a - cutCheckCmdLineFlag et dérivés

Vérifie si le paramètre flag_name a été passé dans les arguments du programme.

CUTBoolean cutCheckCmdLineFlag ( const int argc, const char ** argv, const char * flag_name ) ;

Vérifie si le paramètre flag_name a été passé dans les arguments du programme et met dans val la valeur de cet argument.

CUTBoolean cutGetCmdLineArgumenti   ( const int argc, const char ** argv, const char * arg_name,
int   * val ) ;
CUTBoolean cutGetCmdLineArgumentf   ( const int argc, const char ** argv, const char * arg_name,
float * val ) ;
CUTBoolean cutGetCmdLineArgumentstr ( const int argc, const char ** argv, const char * arg_name,
char ** val ) ;

IX-B-4 - Comparaison

IX-B-4-a - cutComparef et dérivés

Compare les deux tableaux de longueur len.

CUTBoolean cutComparef  ( const      float    * reference, const      float    * data,
const unsigned int len ) ;
CUTBoolean cutComparei  ( const       int     * reference, const       int     * data,
const unsigned int len ) ;
CUTBoolean cutCompareub ( const unsigned char * reference, const unsigned char * data,
const unsigned int len ) ;

IX-B-4-b - cutComparefe

Compare les deux tableaux de longueur len, en autorisant une approximation epsilon.

CUTBoolean cutComparefe ( const float * reference, const float * data, const unsigned int len,
const float epsilon ) ;

IX-B-5 - Gestion du temps

CUTBolean cutCreateTimer   (       unsigned int * name ) ;
CUTBolean cutDeleteTimer   (       unsigned int   name ) ;
CUTBolean cutStartTimer    ( const unsigned int   name ) ;
CUTBolean cutStopTimer     ( const unsigned int   name ) ;
CUTBolean cutResetTimer    ( const unsigned int   name ) ;
CUTBolean cutGetTimerValue ( const unsigned int   name ) ;

Leur utilisation sera explicitée dans l'exemple.

IX-C - Un exemple complet

int main()
{
float *h_A, *h_B, *h_C_GPU;
float *d_A, *d_B, *d_C;
unsigned int hTimer;
int i, j;

int VECTOR_N  = 512;
int    DATA_N = VECTOR_N;
int   DATA_SZ = DATA_N    * sizeof(float);
int RESULT_SZ = VECTOR_N  * sizeof(float);
CUT_DEVICE_INIT();
CUT_SAFE_CALL( cutCreateTimer(& hTimer) );

printf("Initializing data...\n");
printf("...allocating CPU memory.\n");
h_A     = (float *) malloc(DATA_SZ);
h_B     = (float *) malloc(DATA_SZ);
h_C_GPU = (float *) malloc(RESULT_SZ);

printf("...allocating GPU memory.\n");
CUDA_SAFE_CALL( cudaMalloc( (void **) & d_A, DATA_SZ)   );
CUDA_SAFE_CALL( cudaMalloc( (void **) & d_B, DATA_SZ)   );
CUDA_SAFE_CALL( cudaMalloc( (void **) & d_C, RESULT_SZ) );

printf("...generating input data in CPU mem.\n");
srand(153);
for(i = 0 ; i < DATA_N ; i++)
{
h_A[i] = RandFloat(0.0f, 1.0f);
h_B[i] = RandFloat(0.0f, 1.0f);
}

printf("...copying input data to GPU mem.\n");
CUDA_SAFE_CALL( cudaMemcpy(d_A, h_A, DATA_SZ, cudaMemcpyHostToDevice) );
CUDA_SAFE_CALL( cudaMemcpy(d_B, h_B, DATA_SZ, cudaMemcpyHostToDevice) );
printf("Data init done.\n");

printf("Executing GPU kernel...\n");
CUDA_SAFE_CALL( cudaThreadSynchronize() );
CUT_SAFE_CALL( cutResetTimer(hTimer) );
CUT_SAFE_CALL( cutStartTimer(hTimer) );
scalarProdGPU <<<1, VECTOR_N>>> (d_C, d_A, d_B, VECTOR_N);
CUT_CHECK_ERROR("scalarProdGPU() execution failed\n");
CUDA_SAFE_CALL ( cudaThreadSynchronize() );

printf("Reading back GPU result...\n");
CUDA_SAFE_CALL( cudaMemcpy(h_C_GPU, d_C, RESULT_SZ, cudaMemcpyDeviceToHost) );
CUT_SAFE_CALL ( cutStopTimer(hTimer) );
printf("GPU time: %f msecs.\n", cutGetTimerValue(hTimer));

printf("Shutting down...\n");
CUDA_SAFE_CALL( cudaFree(d_C) );
CUDA_SAFE_CALL( cudaFree(d_B)   );
CUDA_SAFE_CALL( cudaFree(d_A)   );
free(h_C_GPU);
free(h_B);
free(h_A);
CUT_SAFE_CALL( cutDeleteTimer(hTimer) );
}

Ceci ne présente que l'utilisation de fonctions et de macros CUTIL. Des fonctions utilisées sont définies dans le fichier completet dans les kernels.

Les fonctions et macros CUTIL sont définies dans le fichier cutil.h.

X - CUDPP

X-A - Qu'est-ce que CUDPP ?

CUDA Data Parallel Primitives.

Il s'agit d'une autre API C. Elle est donc orientée vers les algorithmes de traitement parallèle des données, comme la réduction parallèle et le tri parallèle des listes. Ces algorithmes, primitifs, sont nécessaires à bien d'autres algorithmes, comme la compression de données, sans être forcément très faciles à implémenter efficacement.

CUDPP se base, comme CuFFT et FFTW pour les transformées de Fourier, sur le principe des plans.

La majorité des développeurs utilise un framework (Qt, GTK+), ou une librairie pluripotente (Boost) pour s'aider dans leurs développements. Ces librairies permettent d'effectuer rapidement, efficacement et simplement des pans entiers de nos programmes, comme le tri. CUDPP effectue les mêmes opérations, mais sur un GPU, à la manière de CuBLAS et de CuFFT : l'utilisateur ne voit presque pas la différence entre GPU et CPU.

Cette librairie est aussi disponible pour Java.

X-B - Compilation

CUDPP n'est pas développé par NVIDIA, ni supporté, par conséquent, non distribué avec le SDK de CUDA. Il vous faudra donc le télécharger puis le compiler. Cette opération prend un peu de temps, environ 15 minutes sur une configuration assez musclée.

Cependant, CUDPP dépend de CUTIL, une librairie développée par NVIDIA, mais non supportée, qui est censée aider au développement avec CUDA. Cette librairie est incluse et doit être compilée avant utilisation.

Commencez par récupérer la librairie sur le site officiel : GPGPU CUDPP. La version actuelle est la 1.0 alpha.

Une fois l'archive récupérée, décompressez-la dans un répertoire, peu importe où, tant qu'il n'y a pas d'espace dans le chemin.

La suite dépend de votre environnement. Vous devrez avoir installé CUDA correctement et configuré l'environnement pour son utilisation.

CUDPP est censé compiler avec CUDA 2.0 et être incompatible avec les autres versions. Or, j'ai pu le compiler sans problème avec CUDA 2.2. Si vous avez des problèmes, il vous suffit de revenir à une version plus ancienne de CUDA.

X-B-1 - Visual Studio

Dans le dossier common, ouvrez les solutions cutil.sln et paramgl.sln. Compilez les configurations Debug et Release.

Ensuite, dans le dossier cudpp, ouvrez la solution cudpp.sln. Sélectionnez la configuration qui vous intéresse (un conseil : compilez au moins une version Debug et une version Release, pour faciliter le débogage de votre application), puis lancez la compilation, comme pour toute solution.

Ensuite, il faut installer ces librairies. Copiez les fichiers .lib du dossier lib dans votre dossier CUDA\lib, les .h du dossier inc dans CUDA\include et les fichiers .dll du dossier lib dans votre dossier CUDA\bin.

Si vous n'arrivez pas à compiler correctement ces librairies, téléchargez cette archive, qui les contient toutes (vous n'avez qu'à l'extraire dans le dossier Visual Studio 9.0\VC). Les fichiers compilés sont prévus pour Visual Studio 9 et CUDA 2.2, toute utilisation avec d'autres versions n'est pas garantie de fonctionnement.

L'utilisation des compilateurs Intel est déconseillée dans ce cas, car nvcc n'a pas encore été testé avec lui ! Le résultat de la compilation peut être très aléatoire !

Pendant la compilation, vous verrez beaucoup d'avertissements : ils ne sont que les signes de l'optimiseur du compilateur, qui optimise tout ce qu'il peut, surtout du côté GPU.

X-B-2 - GCC

Dans un shell Linux (sh, par exemple), lancez ces quelques commandes.

( cd common ; make cuda-install=/usr/local/cuda )
( cd cudpp  ; make cuda-install=/usr/local/cuda )

Elles compilent et installent le tout automatiquement, dans le dossier de CUDA.

X-C - La découverte par l'exemple

Cet exemple est celui livré avec la librairie CUDPP, sous le nom de simpleCUDPP.

1. void runTest(int argc, char ** argv)

2. {

3. CUT_DEVICE_INIT();

5. unsigned int numElements = 32768;

6. unsigned int memSize = sizeof( float) * numElements;

8. // allocate host memory

9. float* h_idata = (float*) malloc( memSize);

10. // initalize the memory

11. for (unsigned int i = 0; i < numElements; ++i)

12. {

13. h_idata[i] = (float) (rand() & 0xf);

14. }

16. // allocate device memory

17. float* d_idata;

18. CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata, memSize));

19. // copy host memory to device

20. CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, memSize,

21. cudaMemcpyHostToDevice) );

23. // allocate device memory for result

24. float* d_odata;

25. CUDA_SAFE_CALL( cudaMalloc( (void**) &d_odata, memSize));

27. CUDPPConfiguration config;

28. config.op = CUDPP_ADD;

29. config.datatype = CUDPP_FLOAT;

30. config.algorithm = CUDPP_SCAN;

31. config.options = CUDPP_OPTION_FORWARD | CUDPP_OPTION_EXCLUSIVE;

33. CUDPPHandle scanplan = 0;

34. CUDPPResult result = cudppPlan(&scanplan, config, numElements, 1, 0);

36. if (CUDPP_SUCCESS != result)

37. {

38. printf("Error creating CUDPPPlan\n");

39. exit(-1);

40. }

42. // Run the scan

43. cudppScan(scanplan, d_odata, d_idata, numElements);

45. // allocate mem for the result on host side

46. float* h_odata = (float*) malloc( memSize);

47. // copy result from device to host

48. CUDA_SAFE_CALL( cudaMemcpy( h_odata, d_odata, memSize,

49. cudaMemcpyDeviceToHost) );

50. // compute reference solution

51. float* reference = (float*) malloc( memSize);

52. computeSumScanGold( reference, h_idata, numElements, config);

54. // check result

55. CUTBoolean res = cutComparef( reference, h_odata, numElements);

56. printf( "Test %s\n", (1 == res) ? "PASSED" : "FAILED");

58. result = cudppDestroyPlan(scanplan);

59. if (CUDPP_SUCCESS != result)

60. {

61. printf("Error destroying CUDPPPlan\n");

62. exit(-1);

63. }

65. free( h_idata);

66. free( h_odata);

67. free( reference);

68. CUDA_SAFE_CALL(cudaFree(d_idata));

69. CUDA_SAFE_CALL(cudaFree(d_odata));

70. }

L'entièreté du fichier est disponible ici et les fonctions utilisées non définies dedans le sont dans ce fichier.

X-D - Décryptage de l'exemple

X-D-1 - Initialisation

Tout d'abord, qu'est-il censé effectuer ? Un scan, aussi dénommé somme partielle, somme des préfixes ou réduction des préfixes. Mais qu'est-ce donc ?

Prenez une liste. Appliquez-lui un scan. Il en résultera une autre liste, dont chaque élément est la somme des précédents.

Liste d'origine

1
2
4
8
16
14
2004
0

Liste scannée

1
3
7
15
31
45
2049
2049

1. CUT_DEVICE_INIT();

3. unsigned int numElements = 32768;

4. unsigned int memSize = sizeof(float) * numElements;

7. float* h_idata = (float *) malloc(memSize);

9. for (unsigned int i = 0; i < numElements; ++i)

10. {

11. h_idata[i] = (float) (rand() & 0xf);

12. }

CUTIL initialise un périphérique, on spécifie le nombre d'éléments à traiter puis la mémoire utilisée pour stockerces éléments. On initialise le tableau qui va les contenir et on le remplit de données aléatoires (entre 0 et 15).

1. float* d_idata;

2. CUDA_SAFE_CALL( cudaMalloc( (void **) & d_idata, memSize));

4. CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, memSize,

5. cudaMemcpyHostToDevice) );

8. float* d_odata;

9. CUDA_SAFE_CALL( cudaMalloc( (void **) & d_odata, memSize));

On crée un tableau semblable sur le GPU et on y transfère les données. On prépare aussi un emplacement pour le résultat.

Maintenant, nous pouvons enfin commencer à parler de CUDPP.

X-D-2 - Les plans

Nous devons commencer par configurer CUDPP correctement pour scanner notre tableau. La configuration des algorithmes de CUDPP repose sur le concept des plans.

Un plan est, tout comme dans CUFFT, une structure qui contient les informations sur les algorithmes à employer pour le calcul.

Ici, nous allons effectuer un scan sur des flottants, exclusivement vers l'avant, d'un nombre d'éléments connu d'avance.

1. CUDPPConfiguration config;

2. config.op = CUDPP_ADD;

3. config.datatype = CUDPP_FLOAT;

4. config.algorithm = CUDPP_SCAN;

5. config.options = CUDPP_OPTION_FORWARD | CUDPP_OPTION_EXCLUSIVE;

7. CUDPPHandle scanplan = 0;

8. CUDPPResult result = cudppPlan(&scanplan, config, numElements, 1, 0);

10. if (CUDPP_SUCCESS != result)

11. {

12. printf("Error creating CUDPPPlan\n");

13. exit(-1);

14. }

X-D-3 - Le calcul

1. cudppScan(scanplan, d_odata, d_idata, numElements);

CUDPP effectue le scan, enfin.

1. float* h_odata = (float*) malloc( memSize);

3. CUDA_SAFE_CALL( cudaMemcpy( h_odata, d_odata, memSize,

4. cudaMemcpyDeviceToHost) );

6. float* reference = (float*) malloc( memSize);

7. computeSumScanGold( reference, h_idata, numElements, config);

10. CUTBoolean res = cutComparef( reference, h_odata, numElements);

11. printf( "Test %s\n", (1 == res) ? "PASSED" : "FAILED");

On récupère le résultat du calcul, on effectue le même calcul en local (sur le CPU) puis on compare les résultats : s'ils sont identiques, le test est réussi.

1. result = cudppDestroyPlan(scanplan);

2. if (CUDPP_SUCCESS != result)

3. {

4. printf("Error destroying CUDPPPlan\n");

5. exit(-1);

6. }

8. free(h_idata);

9. free(h_odata);

10. free(reference);

11. CUDA_SAFE_CALL(cudaFree(d_idata));

12. CUDA_SAFE_CALL(cudaFree(d_odata));

13. }

On détruit le plan, puis on libère la mémoire allouée.

Maintenant, vous pouvez voir la simplicité de CUDPP : on alloue de la mémoire, on la remplit, on configure CUDPP en lui disant ce qu'il doit faire, puis il le fait et on récupère le résultat. C'est tout. Et c'est optimisé.

X-E - Gestion des erreurs

Chaque fonction de CUDPP renvoie un CUDPPResult, qui peut prendre quelques valeurs, en fonction de l'erreur qui s'est déroulée.

Code

Signification

CUDPP_SUCCESS

Pas d'erreur.

CUDPP_ERROR_INVALID_HANDLE

CUDPPHandle spécifié invalide.

CUDPP_ERROR_ILLEGAL_CONFIGURATION

La configuration demandée est impossible (illogique ou invalide).

CUDPP_ERROR_UNKNOWN

Erreur inconnue ou impossible à tracer.

X-F - Options de l'algorithme

Tout algorithme de CUDPP doit être préalablement configuré, grâce à une structure de type CUDPPConfiguration. Cette structure est ainsi définie.

Définition de CUDPPConfiguration

struct CUDPPConfiguration
{
CUDPPAlgorithm algorithm;
CUDPPOperator  op;
CUDPPDatatype  datatype;
unsigned int   options;
};

Nous allons étudier les valeurs de chacune de ces variables.

X-F-1 - Algorithme

Premièrement, nous utilisons une variable de type CUDPPAlgorithm pour configurer un algorithme. Il s'agit, justement, de l'algorithme à employer.

- CUDPP_SCAN

- CUDPP_SEGMENTED_SCAN

- CUDPP_COMPACT

- CUDPP_REDUCE

- CUDPP_SORT_RADIX

- CUDPP_SORT_RADIX_GLOBAL

- CUDPP_SPMVMULT

- CUDPP_SORT_INVALID

Évidemment, cet algorithme devra correspondre à la fonction appelée plus tard avec cette configuration.

X-F-2 - Opérateur

Deuxièmement, vous pouvez spécifier un opérateur, l'opération qui sera effectuée, uniquement lors d'un scan. Pour les autres algorithmes, laissez ceci vide.

Tous ces opérateurs sont binairement associables.

- CUDPP_ADD

- CUDPP_MULTIPLY

- CUDPP_MIN

- CUDPP_MAX

X-F-3 - Type

Troisièmement, vous devez définir le type des données en entrée et en sortie. Voici un tableau récapitulatif de ces valeurs, avec leur équivalent en C.

CUDPPDatatype

C

CUDPP_CHAR

char

CUDPP_UCHAR

unsigned char

CUDPP_INT

int

CUDPP_UINT

unsigned int

CUDPP_FLOAT

float

X-F-4 - Configuration

La configuration d'un algorithme se fait, finalement, aussi grâce à un entier. Cet entier est constitué d'une combinaison de paramètres, grâce à l'opérateur pipe (|).

Nom

Action

Valeur

CUDPP_OPTION_FORWARD

L'algorithme se déroulera du début vers la fin des données.

0x1

CUDPP_OPTION_BACKWARD

L'algorithme se déroulera de la fin vers le début des données.

0x2

CUDPP_OPTION_EXCLUSIVE

L'algorithme inclura tous les éléments précédents sauf l'actuel (uniquement pour des scans).

0x4

CUDPP_OPTION_INCLUSIVE

L'algorithme inclura tous les éléments précédents et l'actuel (uniquement pour des scans).

0x8

CUDPP_OPTION_CTA_LOCAL

L'algorithme sera effectué sans communication entre les blocs (ne fonctionne que pour les tris, pour le moment).

0x10

X-G - Fonctions de plans

CUDPPResult cudppPlan(CUDPPHandle        * planHandle,
CUDPPConfiguration config,
size_t             n,
size_t             rows,
size_t             rowPitch);

CUDPPResult cudppDestroyPlan(CUDPPHandle plan);

cudppPlan() crée un plan planHandler, qui respectera la configuration config, qui traitera n éléments répartis en rows colonnes de rowPitch éléments chacune.

Le plan sera détruit par la fonction cudppDesroyPlan() dès qu'il ne sera plus utilisé.

X-H - Fonctions de calcul

X-H-1 - Scan

CUDPPResult cudppScan(CUDPPHandle   planHandle,
void        * d_out,
const void  * d_in,
size_t        numElements);

Les données à traiter se situent sur le GPU, à l'emplacement pointé par d_in ; la sortie se situe aussi sur le GPU, à l'emplacement pointé par d_out. La fonction utilisera les options prévues dans le plan planHandle et scannera numElements éléments.

CUDPPResult cudppSegmentedScan(CUDPPHandle          planHandle,
void               * d_out,
const void         * d_idata,
const unsigned int * d_iflags,
size_t               numElements);

La fonction traitera les segments précisés dans d_iflags de d_idata.

CUDPPResult cudppMultiScan(CUDPPHandle planHandle,
void        *d_out,
const void  *d_in,
size_t      numElements,
size_t      numRows);

La fonction effectuera numRows scans sur les données avant d'écrire le résultat.

CUDPPResult cudppCompact(CUDPPHandle          planHandle,
void               * d_out,
size_t             * d_numValidElements,
const void         * d_in,
const unsigned int * d_isValid,
size_t               numElements);

Cette fonction compacte un tableau, avec les éléments d_isValid définis comme valides, les autres étant à éjecter. d_numValidElements sera le nombre d'éléments effectivement gardés.

Par exemple :

Exemple de compactage

d_in        = [ a b c d e f ]
deviceValid = [ 1 0 1 1 0 1 ]
d_out       = [ a c d f ]

X-H-2 - Tri

CUDPPResult cudppSort(CUDPPHandle planHandle,
void        *d_out,
const void  *d_in,
size_t      numElements);

Trie le tableau.

XI - Thrust

XI-A - Qu'est-ce que Thrust ?

Il s'agit d'une bibliothèque de templates à utiliser avec CUDA. En tant que telle, elle ne nécessite pas de compilation : pour l'installer, téléchargez la dernière version sur son site web et extrayez l'archive dans votre dossier d'includes de CUDA.

Thrust dispose quand même de fonctions : des fonctions inline. Cela permet de gagner un certain temps à l'exécution : un appel de fonction produit un remplissage des registres du CPU puis un saut vers la fonction, qui resautera vers l'emplacement précédent, ce qui prend un certain temps. En contre partie, votre exécutable sera plus lourd (comme si vous aviez inclus une nouvelle librairie en statique).

Thrust se base sur la STL pour fournir ses fonctions : il arrive assez fréquemment que des fonctions aient le même nom. C'est pourquoi toutes les fonctions de Thrust sont dans le namespace thrust.

XI-B - Un premier exemple : la version

#include <thrust/version.h>
#include <iostream>

int main(void)
{
int major = (THRUST_VERSION) / 100000;
int minor = (THRUST_VERSION  / 100) % 1000;
std::cout << "Thrust v" << major << "." << minor << std::endl;
return 0;
}

Ce code montre comment récupérer la version de Thrust, tout simplement.

XI-C - Les vecteurs

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>

int main(void)
{
thrust::host_vector<int> H(4);

H[0] = 14;
H[1] = 20;
H[2] = 38;
H[3] = 46;

std::cout << "H has size " << H.size() << std::endl;

for(int i = 0; i < H.size(); i++)
std::cout << "H[" << i << "] = " << H[i] << std::endl;

H.resize(2);
std::cout << "H now has size " << H.size() << std::endl;

thrust::device_vector<int> D = H;

D[0] = 99;
D[1] = 88;

for(int i = 0; i < D.size(); i++)
std::cout << "D[" << i << "] = " << D[i] << std::endl;

return 0;
}

Le terme vecteur est ici à comprendre au sens de la STL et non des mathématiques : il ne s'agit pas d'un ensemble de nombres, mais bien d'un tableau dynamique.

Il existe deux types de vecteurs : l'un pour l'hôte (host_vector), l'autre pour le périphérique (device_vector).

Les copies de l'hôte vers le périphérique sont simplifiées à l'extrême : on initialise un vecteur de chaque type, et on affecte la valeur de l'un à l'autre. C'est tout.

Chaque modification d'un vecteur sur le périphérique fait appel à cudaMemcpy(), ils sont donc fort coûteux. Il faut les limiter au maximum.

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/copy.h>
#include <thrust/fill.h>
#include <thrust/sequence.h>
#include <iostream>

int main(void)
{
thrust::device_vector<int> D(10, 1);

thrust::fill(D.begin(), D.begin() + 7, 9);

thrust::host_vector<int> H(D.begin(), D.begin() + 5);

thrust::sequence(H.begin(), H.end());

thrust::copy(H.begin(), H.end(), D.begin());

for(int i = 0; i < D.size(); i++)
std::cout << "D[" << i << "] = " << D[i] << std::endl;
return 0;
}

Ici, on initialise un vecteur sur le périphérique, on remplit les 7 premiers éléments avec 9. On crée un autre vecteur, sur l'hôte, qui contient les 5 premières valeurs du premier vecteur. L'hôte est rempli avec des valeurs incrémentant à partir de 1. On copie l'hôte de H.begin() à H.end() sur D.begin()

XI-D - Itérateurs et pointeurs

Vous n'êtes probablement pas sans savoir que les fonctions begin() et end() que nous venons d'utiliser sont des itérateurs : ils reprennent le nom utilisé pour la STL. Ces fonctions sont implémentées comme des templates : une seule et unique implémentation pour gérer l'hôte et le périphérique. Leurs types de renvoi sont quand même différents selon le milieu où elles doivent être exécutées.

Vous pouvez aussi utiliser des pointeurs créés directement par CUDA avec Thrust.

size_t N = 10;

int * raw_ptr;
cudaMalloc( (void **) & raw_ptr, N * sizeof(int));

thrust::device_ptr<int> dev_ptr(raw_ptr);

thrust::fill(dev_ptr, dev_ptr + N, (int) 0);

Vous pouvez aussi récupérer un pointeur CUDA depuis un pointeur Thrust.

size_t N = 10;

thrust::device_ptr<int> dev_ptr = thrust::device_malloc<int>(N);

int * raw_ptr = thrust::raw_pointer_cast(dev_ptr);

Thrust peut aussi être compatible avec les conteneurs de la STL, quand il n'en dispose pas.

#include <thrust/device_vector.h>
#include <thrust/copy.h>
#include <list>
#include <vector>

int main(void)
{
std::list<int> stl_list;
stl_list.push_back(10);
stl_list.push_back(20);
stl_list.push_back(30);
stl_list.push_back(40);

thrust::device_vector<int> D(stl_list.begin(), stl_list.end());

std::vector<int> stl_vector(D.size());
thrust::copy(D.begin(), D.end(), stl_vector.begin());
return 0;
}

XI-E - Algorithmes

Thrust dispose d'un certain nombre d'algorithmes parallèles, dont une certaine majorité possède un équivalent dans la STL. Les deux versions (STL et Thrust) ont alors le même nom.

Chacun des algorithmes de Thrust est implémenté pour l'hôte et pour le périphérique. Suivant le conteneur ou itérateur qui est passé à l'algorithme, la version hôte ou périphérique sera utilisée.

La majorité des algorithmes de Thrust requièrent que les itérateurs soient dans le même milieu (hôte ou périphérique), à quelques exceptions logiques près, comme copy().

XI-E-1 - Transformations

Ces algorithmes appliquent une même opération à tous les éléments d'une entrée et les stockent à un autre endroit.

#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/sequence.h>
#include <thrust/copy.h>
#include <thrust/fill.h>
#include <thrust/replace.h>
#include <thrust/functional.h>
#include <iostream>

int main(void)
{
thrust::device_vector<int> X(10);
thrust::device_vector<int> Y(10);
thrust::device_vector<int> Z(10);

// remplit X par des valeurs de 1 à 10
thrust::sequence(X.begin(), X.end());

// Y = -X
thrust::transform(X.begin(), X.end(), Y.begin(), thrust::negate<int>());

// remplit Z avec des 2
thrust::fill(Z.begin(), Z.end(), 2);

// Y = X mod 2
thrust::transform(X.begin(), X.end(), Z.begin(), Y.begin(), thrust::modulus<int>());

// remplace tous les 1 de Y par  10
thrust::replace(Y.begin(), Y.end(), 1, 10);

// copie le contenu de Y dans le flux de sortie cout
thrust::copy(Y.begin(), Y.end(), std::ostream_iterator<int>(std::cout, "\n"));

return 0;
}

thrust::modulus et thrust::negate sont implémentés comme des foncteurs, dans <thrust/functional.h>.

Nous allons poursuivre avec une implémentation personnelle d'une fonction de BLAS : SAXPY, qui effectue les opérations de la forme .

include <thrust/transform.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/functional.h>
#include <iostream>
#include <iterator>
#include <algorithm>

struct saxpy_functor
{
const float a;
saxpy_functor(float _a) : a(_a) {}
__host__ __device__ float operator() (const float& x, const float& y) const
{
return a * x + y;
}
};

void saxpy_fast(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
{
thrust::transform(X.begin(), X.end(), Y.begin(), Y.begin(), saxpy_functor(A));
}

void saxpy_slow(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
{
thrust::device_vector<float> temp(X.size());
thrust::fill(temp.begin(), temp.end(), A);
thrust::transform(X.begin(), X.end(), temp.begin(), temp.begin(), thrust::multiplies<float>());
thrust::transform(temp.begin(), temp.end(), Y.begin(), Y.begin(), thrust::plus<float>());
}

Vous pouvez voir deux implémentations : une rapide et une lente. La première effectue lectures et écritures. La seconde effectue lectures et écritures.

Vu que SAXPY est limité par la mémoire et non par la puissance de calcul, la première est plus rapide que la seconde.

La technique utilisée s'appelle la réduction des kernels, car, pour obtenir de meilleures performances, dans ce cas, vous devez effectuer le moins possible de lectures et d'écritures en mémoire.

Cette manière de faire (avec thrust::transform) limite les arguments : un ou deux en entrée. Si vous avez besoin de plus d'arguments, voyez thrust::counting_iterator et thrust::for_each tels qu'utilisés dans ce nouvel exemple.

#include <thrust/for_each.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/generate.h>
#include <thrust/iterator/iterator_traits.h>
#include <thrust/device_ptr.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>
#include <cstdlib>

struct complicated_functor
{
const float * A;
const float * B;
const float * C;
float * D;

complicated_functor(const float * _A, const float * _B, const float * _C, float * _D):
A(_A), B(_B), C(_C), D(_D) {}

template <typename IndexType> __host__ __device__
void operator()(const IndexType & i)
{
D[i] = A[i] + B[i] * C[i];
}
};

int main(void)
{
thrust::device_vector<float> A(5);
thrust::device_vector<float> B(5);
thrust::device_vector<float> C(5);
thrust::device_vector<float> D(5);

A[0] = 3;  B[0] = 6;  C[0] = 2;
A[1] = 4;  B[1] = 7;  C[1] = 5;
A[2] = 0;  B[2] = 2;  C[2] = 7;
A[3] = 8;  B[3] = 1;  C[3] = 4;
A[4] = 2;  B[4] = 8;  C[4] = 3;

thrust::experimental::counting_iterator<int, thrust::random_access_device_iterator_tag> begin(0);

complicated_functor op(thrust::raw_pointer_cast(&A[0]),
thrust::raw_pointer_cast(&B[0]),
thrust::raw_pointer_cast(&C[0]),
thrust::raw_pointer_cast(&D[0]));

thrust::for_each(begin, begin + 5, op);

for(int i = 0; i < 5; i++)
std::cout << A[i] << " + " << B[i] << " * " << C[i] << " = " << D[i] << std::endl;
}

XI-E-2 - Réductions

Un algorithme de réduction utilise une opération binaire pour réduire une entrée à une valeur. Voici la manière d'obtenir la somme d'un tableau avec Thrust.

int sum = thrust::reduce(D.begin(), D.end(), (int) 0, thrust::plus<int>());

Cette manière de procéder convient généralement, mais pas toujours, par exemple pour compter. Thrust, comme la STL, fournit quelques autres fonctions, pour votre facilité.

#include <thrust/count.h>
#include <thrust/device_vector.h>

thrust::device_vector<int> vec(5,0);
vec[1] = 1;
vec[3] = 1;
vec[4] = 1;

int result = thrust::count(vec.begin(), vec.end(), 1);
// result == 3

On peut utiliser ces réductions pour calculer la norme d'un vecteur.

#include <thrust/transform_reduce.h>
#include <thrust/functional.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <cmath>

template <typename T>
struct square
{
__host__ __device__
T operator()(const T& x) const {
return x * x;
}
};

int main(void)
{
float x[4] = {1.0, 2.0, 3.0, 4.0};

thrust::device_vector<float> d_x(x, x + 4);

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

float norm = std::sqrt( thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op) );

std::cout << norm << std::endl;

return 0;
}

En fusionnant la multiplication avec l'addition, nous avons un kernel très optimisé, autant qu'un kernel écrit à la main et patiemment optimisé.

XI-E-3 - Scans

#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::inclusive_scan(data, data + 6, data);
// data == {1, 1, 3, 5, 6, 9}

#include <thrust/scan.h>
int data[6] = {1, 0, 2, 2, 1, 3};
thrust::exclusive_scan(data, data + 6, data);
// data == {0, 1, 1, 3, 5, 6}

XI-E-4 - Tri

Thrust, une fois de plus, se calque sur la STL pour fournir ses algorithmes de tri.

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <cstdlib>

int main(void)
{
thrust::host_vector<int> h_vec(20);
thrust::generate(h_vec.begin(), h_vec.end(), rand);

thrust::device_vector<int> d_vec = h_vec;

thrust::sort(d_vec.begin(), d_vec.end());

return 0;
}

De plus, Thrust propose le tri par clé.

const int N = 6;
int    keys[N] = {  1,   4,   2,   8,   5,   7};
char values[N] = {'a', 'b', 'c', 'd', 'e', 'f'};
thrust::sort_by_key(keys, keys + N, values);
// keys   == {  1,   2,   4,   5,   7,   8}
// values == {'a', 'c', 'b', 'e', 'f', 'd'}

Et selon un critère défini par l'utilisateur.

const int N = 6;
int A[N] = {1, 4, 2, 8, 5, 7};
thrust::stable_sort(A, A + N, thrust::greater<int>());
// A == {8, 7, 5, 4, 2, 1}

XII - Fonctions mathématiques

Ces fonctions sont des équivalents des fonctions proposées par la librairie standard du C. Cependant, toutes ne le sont pas. Si ces fonctions sont utilisées dans du code pour l'hôte, alors une implémentation CUDA est utilisée. Sinon, c'est l'implémentation de la librairie standard qui est utilisée.

Toutes ces fonctions existent en deux versions : simple et double précision. Si la version double précision est appelée avec un argument de simple précision, c'est la version de simple précision qui sera appelée.

Voici les 68 fonctions supportées :

- x+y

- x*y

- x/y

- 1/x

- 1/sqrtf(x)

- rsqrtf(x)

- sqrtf(x)

- cbrtf(x)

- hypotf(x,y)

- expf(x)

- exp2f(x)

- exp10f(x)

- expm1f(x)

- logf(x)

- log2f(x)

- log10f(x)

- log1pf(x)

- sinf(x)

- cosf(x)

- tanf(x)

- sincosf(x,sptr,cptr)

- asinf(x)

- acosf(x)

- atanf(x)

- atan2f(y,x)

- sinhf(x)

- coshf(x)

- tanhf(x)

- asinhf(x)

- acoshf(x)

- atanhf(x)

- powf(x,y)

- erff(x)

- erfcf(x)

- lgammaf(x)

- tgammaf(x)

- fmaf(x,y,z)

- frexpf(x,exp)

- ldexpf(x,exp)

- scalbnf(x,n)

- scalblnf(x,l)

- logbf(x)

- ilogbf(x)

- fmodf(x,y)

- remainderf(x,y)

- remquof(x,y,iptr)

- modff(x,iptr)

- fdimf(x,y)

- truncf(x)

- roundf(x)

- rintf(x)

- nearbyintf(x)

- ceilf(x)

- floorf(x)

- lrintf(x)

- lroundf(x)

- llrintf(x)

- llroundf(x)

- signbit(x)

- isinf(x)

- isnan(x)

- isfinite(x)

- copysignf(x,y)

- fminf(x,y)

- fmaxf(x,y)

- fabsf(x)

- nanf(cptr)

- nextafterf(x,y)

Toutes ces opérations ne sont pas aussi précises l'une que l'autre : la table des précisions est disponible dans le Programming Guide, dans l'annexe B.

XIII - Fonctions atomiques

Ces fonctions ne sont pas disponibles sur tous les périphériques : seulement ceux de Compute Capability supérieure ou égale à 1.1. Ceci se résume ainsi : les cartes GeForce 9000 et supérieures, les Quadro à l'exception des FX 4600 et FX 5600 et les Quadro Plex à l'exception des 1000 et 2100 S4.

Qu'est-ce qu'une fonction atomique ?
Une fonction atomique est une opération ou un ensemble d'opérations qui s'exécute entièrement sans que le thread qui l'a demandée ne cède sa place pendant son traitement.
Autrement dit, lors d'un appel à une fonction atomique, aucun autre thread CPU ne peut accéder au même processeur de flux.

Les fonctions atomiques de CUDA se divisent en deux groupes : les fonctions arithmétiques et les fonctions d'opération sur les bits.

Ces fonctions ne peuvent se dérouler qu'en mémoire partagée sur des mots de 32 bits. Seuls les périphériques de Compute Capability de 1.2 ou plus peuvent les effectuer sur des mots de 64 bits. Les périphériques qui en sont capables sont les GeForce GTX 200, les Tesla C1000 et les Quadro FX 5800. Ces périphériques sont aussi les seuls à pouvoir effectuer des calculs à double précision (FP64).

XIII-A - Arithmétiques

int atomicAdd(int * address, int val);
unsigned int atomicAdd(unsigned int * address, unsigned int val);
unsigned long long int atomicAdd(unsigned long long int * address, unsigned long long int val);

Cette fonction lit la variable stockée à l'adresse address et lui ajoute val. La fonction renvoie l'ancienne valeur, avant comparaison.

int atomicSub(int * address, int val);
unsigned int atomicSub(unsigned int * address, unsigned int val);

Cette fonction lit la variable stockée à l'adresse address et lui retire val. La fonction renvoie l'ancienne valeur, avant comparaison.

int atomicExch(int * address, int val);
unsigned int atomicExch(unsigned int * address, unsigned int val);
unsigned long long int atomicExch(unsigned long long int * address, unsigned long long int val);
float atomicExch(float * address, float val);

Cette fonction lit la variable stockée à l'adresse address et lui met val comme valeur. La fonction renvoie l'ancienne valeur, avant remplacement.

int atomicMin(int* address, int val);
unsigned int atomicMin(unsigned int* address, unsigned int val);

Cette fonction lit la variable stockée à l'adresse address et la remplace par le minimum de cette variable et de val. La fonction renvoie l'ancienne valeur, avant calcul.

int atomicMax(int* address, int val);
unsigned int atomicMax(unsigned int* address, unsigned int val);

Cette fonction lit la variable stockée à l'adresse address et la remplace par le maximum de cette variable et de val. La fonction renvoie l'ancienne valeur, avant calcul.

unsigned int atomicInc(unsigned int* address, unsigned int val);

Cette fonction lit la variable stockée à l'adresse address et l'incrémente si elle est supérieure ou égale à val, sinon l'annule. La fonction renvoie l'ancienne valeur, avant calcul.

unsigned int atomicDec(unsigned int* address, unsigned int val);

Cette fonction lit la variable stockée à l'adresse address et la décrémente si elle est supérieure ou égale à val, sinon l'annule. La fonction renvoie l'ancienne valeur, avant calcul.

typedef unsigned int uint;
typedef unsigned long long int ullint,
int atomicCAS(int * address, int compare, int val);
uint atomicCAS(uint * address, uint compare, uint val);
ullint atomicCAS(ullint * address, ullint compare, ullint val);

Cette fonction lit la variable stockée à l'adresse address et la remplace par val si elle est égale à compare. La fonction renvoie l'ancienne valeur, avant calcul.

XIII-B - Bits

Il n'existe pas une aussi grande quantité d'opérations sur les bits. Seules trois sont supportées par CUDA : et (and), ou (or) et ou exclusif (xor). Le premier paramètre de l'opération est la variable dont un pointeur est passé, le second est la variable passée en argument. La fonction renvoie, comme toujours, la valeur sur le GPU avant calcul.

int atomicAnd(int * address, int val);
unsigned int atomicAnd(unsigned int * address, unsigned int val);

int atomicOr(int * address, int val);
unsigned int atomicOr(unsigned int * address, unsigned int val);

int atomicXor(int * address, int val);
unsigned int atomicXor(unsigned int * address, unsigned int val);

XIV - Fonction de synchronisation

Dans le runtime que nous étudions, il n'existe qu'une seule fonction de synchronisation des threads : void __syncthreads();.

Elle synchronise tous les threads d'un bloc. Dès qu'ils sont tous arrivés à ce point dans les calculs, ils recommencent.

Cette synchronisation est principalement utilisée pour coordonner la communication entre les threads d'un même bloc. Quand quelques threads accèdent au même emplacement mémoire, il peut y avoir quelques problèmes, vu que l'on ne peut jamais savoir si un thread a déjà écrit ou lu des données. Ceci peut être évité avec les threads : on peut être sûr que tous les threads sont arrivés à cet endroit, qu'ils ont tous lu ou écrit ce qu'ils devaient.

Cette opération est aussi fort coûteuse en temps : elle prend 4 cycles, si elle ne doit attendre aucun thread. Mais elle peut éviter des erreurs d'exécution, qui doivent être évitées à n'importe quel prix.

XV - Interopérabilité avec DirectX

Des ressources Direct3D peuvent être utilisées avec CUDA, pour que CUDA puisse lire ces données pour les traiter, ou bien pour que CUDA puisse écrire des données qui seront utilisées par Direct3D.

CUDA peut fonctionner avec les versions 9.0 et 10.0 de Direct3D. Pour la première, les fonctions sont préfixées par cudaD3D9. Pour la seconde, par cudaD3D10. C'est le seul élément qui change entre les deux fonctions (à part le fait qu'elles ne sont pas utilisées avec les mêmes versions de D3D). C'est pourquoi je ne préciserai que les fonctions pour D3D9.

Cependant, tout ne peut pas fonctionner comme vous le désirez, mais ces restrictions restent très logiques. Un contexte CUDA ne peut communiquer qu'avec un seul périphérique D3D à la fois et ils doivent correspondre au même GPU. De plus, le périphérique D3D doit être créé avec le drapeau D3DCREATE_HARDWARE_VERTEXPROCESSING.

if( NULL == ( g_pD3D = Direct3DCreate9( D3D_SDK_VERSION ) ) )
return E_FAIL;

D3DPRESENT_PARAMETERS d3dpp;
ZeroMemory( & d3dpp, sizeof(d3dpp) );
d3dpp.Windowed              = TRUE;
d3dpp.BackBufferCount       = 1;
d3dpp.SwapEffect            = D3DSWAPEFFECT_COPY;
d3dpp.hDeviceWindow         = g_handles.hwndText;
d3dpp.BackBufferWidth       = IMAGE_SIZE_X;
d3dpp.BackBufferHeight      = IMAGE_SIZE_Y;
d3dpp.BackBufferFormat      = D3DFMT_UNKNOWN;
d3dpp.PresentationInterval  = D3DPRESENT_INTERVAL_IMMEDIATE;

hr = g_pD3D->CreateDevice (adapter, D3DDEVTYPE_HAL, g_handles.hwndCGH,
D3DCREATE_HARDWARE_VERTEXPROCESSING,
&d3dpp, &g_pD3DDevice);
if( FAILED( hr ) )
return hr;

cudaD3D9SetDirect3DDevice(g_pD3DDevice);

La fonction cudaD3D9SetDirect3DDevice() doit être appelée avant toute autre fonction de CUDA !

Pour associer une ressource DirectX à CUDA, cudaD3D9RegisterResource() est utilisé.

IDirect3DVertexBuffer9 * vbuffer;
cudaD3D9RegisterResource(vbuffer, cudaD3D9RegisterFlagsNone);

IDirect3DSurface9 * surface;
cudaD3D9RegisterResource(surface, cudaD3D9RegisterFlagsNone);

ID3D9Buffer * buffer;
cudaD3D9RegisterResource(buffer, cudaD3D9RegisterFlagsNone);

ID3D9Texture2D* tex2D;
cudaD3D9RegisterResource(tex2D, cudaD3D9RegisterFlagsNone);

Vous pouvez désenregistrer les ressources grâce à cudaD3D9UnregisterVertexBuffer()

cudaD3D9UnregisterVertexBuffer(vbuffer);
cudaD3D9UnregisterVertexBuffer(surface);
cudaD3D9UnregisterVertexBuffer(buffer);
cudaD3D9UnregisterVertexBuffer(tex2D);

Dès qu'une ressource est enregistrée, vous pouvez la lier avec cudaD3D9MapResources() et la délier avec cudaD3D9UnmapResources(). Dès qu'elle est liée, vous pouvez l'utiliser dans un kernel. Pour ce faire, vous aurez besoin de son adresse, retournée par cudaD3D9ResourceGetMappedPointer(). Si vous accédez à une ressource liée par D3D, les résultats sont indéfinis.

ID3D9Buffer * buffer;
cudaD3D9RegisterResource(buffer, cudaD3D9RegisterFlagsNone);

void * devPtr;
cudaD3D9ResourceGetMappedPointer(& devPtr, buffer, 0, 0);

size_t size;
cudaD3D9ResourceGetMappedSize(& size, buffer, 0, 0);

cudaMemset(devPtr, 0, size);

Dans ce dernier exemple, chaque thread accède à un pixel de la surface 2D.

void* devPtr;
cudaD3D9ResourceGetMappedPointer(& devPtr, surface, 0, 0);
size_t pitch;
cudaD3D9ResourceGetMappedPitch(& pitch, 0, surface, 0, 0);
dim3 Db = dim3(16, 16);
dim3 Dg = dim3( (width+Db.x-1) / Db.x, (height+Db.y-1) / Db.y);
mykernel <<< Dg, Db >>> ( (unsigned char *) devPtr, width, height, pitch);

__global__ void mykernel(unsigned char * surface,
int width, int height, size_t pitch)
{
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x >= width || y >= height) return;
float * pixel = (float *) (surface + y * pitch) + 4 * x;
}

XVI - Interopérabilité avec OpenGL

Des ressources OpenGL peuvent être utilisées avec CUDA, pour que CUDA puisse lire ces données pour les traiter, ou bien pour que CUDA puisse écrire des données qui seront utilisées par OpenGL.

Tout comme Direct3D, la fonction cudaGLSetGLDevice() doit être appelée avant toutes les autres fonctions CUDA.

Un buffer peut être enregistré avec la fonction cudaGLRegisterBufferObject() et lié avec cudaGLMapBufferObject().

Vous pouvez délier avec cudaGLUnmapBufferObject() et désenregistrer avec cudaGLUnregisterBufferObject()

GLuint bufferObj;

cudaGLRegisterBufferObject(bufferObj);

float * devPtr;

cudaGLMapBufferObject( (void **) & devPtr, bufferObj);

XVII - Gestion des erreurs

Aucune fonction de CUDA ne renvoie directement ses erreurs. Il faut faire appel à une fonction spéciale, cudaGetLastError(), qui renvoie le code d'erreur, sous forme de cudaError_t.

Toutes les fonctions du runtime renvoient donc un code d'erreur, mais dans le cas d'appels asynchrones, comme le kernel retourne immédiatement, on ne peut pas savoir si une erreur a eu lieu. Vous devez donc utiliser une fonction de synchronisation afin de pouvoir récupérer l'erreur.

Chaque erreur est effacée par la suivante et une réussite renvoie aussi un code d'erreur. Récupérer une erreur remet le compteur à la réussite.

Un kernel ne peut pas retourner de code d'erreur, on doit donc vérifier s'il peut s'exécuter sans problème avant de le lancer.

Il existe une seconde fonction pour gérer les erreurs : cudaGetErrorString(), qui prend en paramètre un code d'erreur, tel que retourné par cudaGetLastError(). Elle retourne, sous la forme d'un const char *, la description de l'erreur.

Code d'erreur

Signification

cudaSuccess

Pas d'erreur.

cudaErrorMissingConfiguration

Configuration manquante.

cudaErrorMemoryAllocation

Erreur lors de l'allocation de la mémoire.

cudaErrorInitializationError

Erreur lors de l'initialisation.

cudaErrorLaunchFailure

Erreur au lancement.

cudaErrorPriorLaunchFailure

Erreur lors d'un lancement précédent.

cudaErrorLaunchTimeout

Temps d'exécution trop long.

cudaErrorLaunchOutOfResources

Plus de ressources disponibles pour l'exécution.

cudaErrorInvalidDeviceFunction

Fonction invalide pour le périphérique.

cudaErrorInvalidConfiguration

Configuration invalide.

cudaErrorInvalidDevice

Périphérique invalide.

cudaErrorInvalidValue

Valeur invalide.

cudaErrorInvalidPitchValue

Valeur de pas invalide.

cudaErrorInvalidSymbol

Symbole invalide.

cudaErrorMapBufferObjectFailed

Erreur lors de l'accès au buffer.

cudaErrorUnmapBufferObjectFailed

Erreur lors de l'accès au buffer.

cudaErrorInvalidHostPointer

Pointeur sur l'hôte invalide.

cudaErrorInvalidDevicePointer

Pointeur sur le périphérique invalide.

cudaErrorInvalidTexture

Texture invalide.

cudaErrorInvalidTextureBinding

Binding de texture invalide.

cudaErrorInvalidChannelDescriptor

Descripteur de chaîne invalide.

cudaErrorInvalidMemcpyDirection

Direction invalide pour memcpy.

cudaErrorAddressOfConstant

Erreur dans l'adresse d'une constante.

cudaErrorTextureFetchFailed

Erreur lors de l'accès à une texture.

cudaErrorTextureNotBound

Texture non liée.

cudaErrorSynchronizationError

Erreur de synchronisation.

cudaErrorInvalidFilterSetting

Paramètre de filtre invalide.

cudaErrorInvalidNormSetting

Paramètre de norme invalide.

cudaErrorMixedDeviceExecution

Exécution mixte sur le périphérique.

cudaErrorCudartUnloading

CUDA runtime en cours de déchargement.

cudaErrorUnknown

Erreur inconnue.

cudaErrorNotYetImplemented

Fonction pas encore implémentée.

cudaErrorMemoryValueTooLarge

Valeur mémoire trop grande.

cudaErrorInvalidResourceHandle

Handle de ressource invalide.

cudaErrorNotReady

Pas encore prêt.

cudaErrorInsufficientDriver

CUDA runtime plus récent que le driver.

cudaErrorSetOnActiveProcess

Set on active process error.

cudaErrorNoDevice

Pas de périphérique disponible.

cudaErrorStartupFailure

Échec du lancement.

cudaErrorApiFailureBase

Erreur dans l'API.

XVIII - Gestion des périphériques

cudaGetDeviceCount() et cudaGetDeviceProperties() permettent d'énumérer les périphériques, les cartes graphiques et leurs capacités.

int deviceCount;
cudaGetDeviceCount(& deviceCount);
int device;

for (device = 0 ; device < deviceCount ; ++device)
{
cudaDeviceProp deviceProp;
cudaGetDeviceProperties(&deviceProp, device);
}

Seuls les périphériques avec un indice de capacité de calcul supérieur ou égal à 1.0 sont retournés, les autres ne pouvant pas être exploités par CUDA.

cudaSetDevice() sert à sélectionner le périphérique qui exécutera les calculs du thread CPU.

cudaSetDevice(device);

Un périphérique doit être choisi avant d'appeler le moindre kernel, soit explicitement avec cudaSetDevice(), soit implicitement (le périphérique 0 est alors choisi). Tout appel ultérieur à cudaSetDevice() échouera.

La fonction cudaGetDeviceProperties() renvoie une structure ainsi définie.

struct cudaDeviceProp
{
char name[256];
size_t totalGlobalMem;
size_t sharedMemPerBlock;
int regsPerBlock;
int warpSize;
size_t memPitch;
int maxThreadsPerBlock;
int maxThreadsDim[3];
int maxGridSize[3];
size_t totalConstMem;
int major;
int minor;
int clockRate;
size_t textureAlignment;
int deviceOverlap;
int multiProcessorCount;
int kernelExecTimeoutEnabled;
}

- name : chaîne ASCII identifiant le périphérique ;

- totalGlobalMem : mémoire globale totale disponible sur le périphérique en bytes ;

- sharedMemPerBlock : mémoire partagée disponible pour un bloc en bytes (cette mémoire est partagée par tous les blocs de threads résidant simultanément sur un multiprocesseur) ;

- regsPerBlock : nombre de registres 32 bits disponibles à un bloc de threads ;

- warpSize : taille d'un warp en threads ;

- memPitch : pas maximal en bytes permis par la copie de mémoire qui implique des régions de mémoire allouées par cudaMallocPitch(); ;

- maxThreadsPerBlock : nombre maximal de threads par bloc ;

- maxThreadsDim[3] : dimensions maximales d'un bloc ;

- maxGridSize[3] : dimensions maximales d'une grille ;

- totalConstMem : total de mémoire constante disponible sur le périphérique ;

- major : numéro de révision majeur de la capacité de calcul ;

- minor : numéro de révision mineur de la capacité de calcul ;

- clockRate : fréquence de l'horloge en kilohertz ;

- textureAlignment : alignement des textures requis, les adresses de base des textures qui sont alignées sur ceci n'ont pas besoin d'appliquer un offset lors des recherches ;

- deviceOverlap : 1 si le périphérique peut copier de la mémoire entre l'hôte et le périphérique en même temps qu'exécuter un kernel, 0 sinon ;

- multiProcessorCount : nombre de multiprocesseurs sur le périphérique ;

- kernelExecTimeoutEnabled : 1 s'il y a une limite de temps pour les kernels, 0 sinon.

Vous pouvez définir vous-même une structure de ce genre et demander à CUDA de choisir le périphérique qui y correspond le mieux, grâce à la fonction cudaChooseDevice( int * dev, const struct cudaDeviceProp * prop ) (elle retourne le périphérique dans dev).

XIX - Conclusions

CUDA n'est qu'une extension du C, à la base. Cette extension ne fait que permettre la communication avec les périphériques.

Tout autour, NVIDIA a bâti un immense empire fonctionnel : cet empire vous facilite fortement la vie, avec quelques types, quelques variables, des lacs entiers de fonctions de copie. Ensuite, il vous offre moult autres outils : les textures, les flux..., qui vous permettent d'être de plus en plus efficaces.

Ensuite viennent les sbires de l'empire : CuBLAS, CuFFT et CUTIL, qui ne sont là que pour vous simplifier encore la vie, sans oublier qu'ils sont très optimisés ! Cela ne gâche rien du plaisir.

Puis les mercenaires, ces étrangers qui ont su se faire une place de choix dans cet immense empire. Il s'agit de CUDPP et de Thrust. Ils sont accompagnés par bien d'autres, tout aussi illustres, mais il faut faire un choix et c'est ceux-là que j'ai retenu.

Finalement, les rats : ces petits outils, tout petits, mais d'une si grande utilité. Ils accélèrent, ils réduisent les erreurs possibles. En un mot : indispensables.

L'univers de CUDA est très riche, très étendu, avec pour objectifs la simplicité et la rapidité. Tous les moyens sont bons pour y arriver.

NVIDIA a aussi eu l'intelligence de permettre à d'autres d'utiliser cette puissance pour encore l'augmenter. Ainsi, CUDA restera encore un certain temps sur le devant de la scène. L'avenir nous dira si OpenCL, un standard proposé par Apple, pourra faire sa place, ou devra s'effacer devant la communauté créée autour de CUDA.

XIX-A - Téléchargements

- Règles de compilation

- Exemple d'utilisation de CuBLAS

- Exemple d'utilisation de CuFFT (une dimension)

- Exemple d'utilisation de CuFFT (deux dimensions)

- Exemple d'utilisation de CUTIL

- kernel utilisé dans l'exemple

- CUDPP compilé pour Visual Studio 2008 (32 et 64 bits)

- Exemple d'utilisation de CUDPP

- Fonctions utilisées dans l'exemple

- Les exemples de Thrust

XIX-B - Voir aussi

- La FAQ GPGPU

- Le forum de Developpez pour vos questions sur CUDA

- La communauté de GPGPU NVIDIA

- Le site web de CUDA

XIX-C - Remerciements

Tout cet article repose sur le Programming Guide diffusé avec le SDK de CUDA, dans ses versions 2.1 et 2.2. Les illustrations (hors formules) proviennent de cet ouvrage. Les formules ont été affichées avec .

Un tout grand merci à raptor70, gorgonite, ange_blond, TanEk et surtout IrmatDen, Matthieu Brucher, ainsi qu'à Étienne Bougoin et Sébastien Flochlay, pour leurs commentaires, encouragements et précisions, sans lesquels l'article ne serait pas celui qu'il est maintenant !

À leurs côtés furent Bovino, mlny84, mais aussi nicorama, Pierre Fauconnier, RideKick, koopajah, jacques_jean, attentifs et courageux relecteurs !