3. Le Modèle Logique

3.1. Grille et bloc

L'essentiel des performances et de l'évolutivité de CUDA réside dans le modèle logique d'exécution qui partitionne les calculs en blocs de threads et les répartit sur l'architecture hautement parallèle du GPU.

La configuration matérielle CUDA joue un rôle important sur les performances d'exécution des kernels comme nous le verrons dans le chapitre suivant.

L'écriture de fonctions parallélisables ne revêt pas de caractère complexe, par contre la répartition des threads qui exécutent les kernels en parallèle pose quelques difficultés lorsque l'on commence à apprendre CUDA.

Lorsqu'un kernel est invoqué, il s'exécute sous forme d'une grille (grid) de blocs de threads parallèles. Dans le cas de l'exemple du chapitre 2, Chaque thread exécutera la partie de code z[i] = x[i] + y[i].

Attention : en prévision de ce que nous verrons dans le chapitre suivant, le nombre maximum de threads dans un bloc est fixé par les caractéristiques de la carte (Compute Capability). Ce nombre maximum de threads sera de 512 ou 1024.

Exemple de grille 2D $(y=2 × x=3)$ avec blocs 2D $(y=3 × x=4)$

Certaines notions sont interchangeables :

• une grille est une matrice de blocs, on peut alors soit parler de grille ou du nombre de blocs de la grille
• un bloc est une matrice (grille) de threads : on peut parler de taille d'un bloc ou du nombre de threads dans un bloc.

Dans l'exemple précédent, on dispose d'une grille de (2 x 3) blocs. Chaque bloc se compose de (3 x 4) threads. Au final on a donc $(2 × 3) × (3 × 4) = 72$ threads à exécuter.

Chaque thread à l'intérieur d'un bloc est identifié par une variable threadIdx (thread index) qui correspond à une structure (dim3) composée de trois champs x, y, z :

typedef struct {
	int x, y, z;
} dim3;

// variable locale au thread au sein du kernel
dim3 threadIdx;

Cette structure matérialise un bloc de threads qui peut être modélisé par une, deux ou trois dimensions. On peut en fonction des besoins n'utiliser :

• qu'une dimension : le thread sera identifié dans le bloc par la coordonnée x (threadIdx.x)
• deux dimensions : le thread sera identifié par x et y (threadIdx.y)
• trois dimensions : le thread sera identifié par x, y et z (threadIdx.z)

Il faudra alors établir une formule de mise en correspondance avec les indices des données en entrée (qui sont des tableaux), le tout en fonction des coordonnées x, y et z.

3.1.1. Exemple 1

Par exemple, si les données en entrée sont codées sous forme d'un tableau à une dimension et que le kernel est exécuté par une grille d'un seul bloc :

• pour un bloc de threads d'une seule dimension ($\,D_x\,$) on utilisera la coordonnée $x$ qui correspond à l'indice i dans le tableau
• pour un bloc de threads à deux dimensions $(\,D_x,\, D_y\,)$ l'identifiant de thread de coordonnées $(x,y)$ sera donné par la formule : $$(\,x + y × D_x\,)$$
• pour un bloc de threads à deux dimensions $(\,D_x,\,D_y,\,D_z\,)$ l'identifiant de thread de coordonnées $(\,x,y,z\,)$ sera donné par la formule : $$(\,x + y × D_x + z × D_x × D_y\,)$$

3.1.2. Exemple 2

Voici un autre exemple simple qui s'appuie sur les données à traiter qui sont des tableaux à deux dimensions de float. Ici on n'a pas à établir de formule de correspondance car on établit une relation directe entre les indices de ligne et de colonne des tableaux avec respectivement y et x :

Remarque, on aurait pu également utiliser des entiers lors de l'appel:

MatAdd <<< 1, N*N >>>(...)

Cependant, ici, on aura pas un bloc de dimension 2D (N,N) mais un bloc 1D de taille N*N ce qui ne permet pas d'établir la correspondance directe avec i et j. Il faudra alors utiliser la formule:

int y = blockIdx.x;
int x = threadIdx.x;

3.2. Appel d'un kernel

L'appel du kernel est réalisé grâce à l'utilisation de l'opérateur d'exécution de configuration :

kernel <<< ... >>>

Cet opérateur peut prendre entre deux et quatre paramètres :

kernel<<< grid, block, shared_mem*, stream* >>>
      (* optionnel)

• grid ( dimensions de la grille) : représente soit un entier, soit un type dim3 et permet de spécifier la dimension de la grille telle que le produit grid.x * grid.y * grid.z correspond au nombre de blocs qui seront créés. Généralement grid.z = 1
• block (dimensions du bloc) : représente soit un entier, soit un type dim3 et définit la taille d'un bloc tel que le produit block.x * block.y * block.z sera égal au nombre de threads par bloc.
• shared_mem (mémoire partagée) : entier naturel qui spécifie la taille en octets de la mémoire partagée allouée par bloc. En général égale à 0 si on a pas de traitement spécifique qui nécessite de la mémoire partagée.
• stream : de type cudaStream_t argument optionnel égal à 0 par défaut qui permet de spécifier un Stream.

Voici quelques exemples d'appel:

kernel<<< 10, 32 >>>(..);
kernel<<< dim3(10), 32 >>>(..);
kernel<<< dim3(10,1,1), dim3(32) >>>(..);

dim3 grid(10);
dim3 block(32);
kernel<<< grid, block >>>(..);

L'appel à une fonction échouera si les paramètres sont erronés ou supérieurs aux tailles maximum définies par la Compute Capability du GPU. Il en va de même si on réserve plus de mémoire partagée que disponible.

Pour traiter, les erreurs, les fonctions CUDA retournent toutes un code d'erreur sauf l'appel au kernel puisqu'un kernel est de type void :

cudaError_t cudaGetLastError(void);
char* cudaGetErrorString(cudaError_t code);

Afin de gérer les erreurs éventuelles qui peuvent survenir il est recommandé de :

• vérifier que l'appel à chaque fonction ne produit pas d'erreur (allocation mémoire, copie des données)
• vérifier que durant l'exécution du kernel aucune erreur n'a été générée (ex: dépassement d'indice, débordement mémoire)

Pour cela on peut définir des macro-instructions en langage C. Le code suivant est issu du package CUME (CUda Made Easy)

3.2.1. Notion de stream

Un Stream représente une série d'opérations à exécuter.

Le stream "0" est complétement synchrone entre host et device et agit comme si on avait inséré l'instruction cudaDeviceSynchronize() avant chaque opération à exécuter sauf si on utilise une instruction asynchrone du type cudaMemcpy*Async

3.2.1.a  Exemple de traitement synchrone

double *cmem = new double [ size ];
cudaMalloc(&gmem, size * sizeof(double) );

// completely synchronous
cudaMemcpy(gmem, cmem, SIZE * sizeof(double), cudaMemcpyHostToDevice);
kernel1<<<grid, block>>>(gmem, ...);
kernel2<<<grid, block>>>(gmem, ...);
cudaMemcpy(cmem, gmem, SIZE * sizeof(double), cudaMemcpyDeviceToHost);

3.2.1.b  Exemples de traitement Asynchrones

Un premier cas de traitement asynchrone apparaît lorsque l'on exécute du code sur le CPU à l'intérieur d'une série d'instructions CUDA :

double *cmem = new double [ size ];
cudaMalloc(&gmem, size * sizeof(double) );

// completely synchronous
cudaMemcpy(gmem, cmem, SIZE * sizeof(double), cudaMemcpyHostToDevice);
kernel1<<<grid, block>>>(gmem, ...);
some_CPU_code();
kernel2<<<grid, block>>>(gmem, ...);
cudaMemcpy(cmem, gmem, SIZE * sizeof(double), cudaMemcpyDeviceToHost);

Un second cas apparaît avec l'utilisation des streams :

cudaStream_t stream1, stream2, ...;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);

cudaMalloc( &dev1, size ) ;
cudaMallocHost( &host1, size ) ; // pinned memory required on host
...
cudaMemcpyAsync( dev1, host1, size, H2D, stream1 ) ;
kernel2 <<< grid, block, 0, stream2 >>> ( dev2, ...) ;
kernel3 <<< grid, block, 0, stream3 >>> ( dev3, ...) ;
cudaMemcpyAsync( host4, dev4, size, D2H, stream4 ) ;
some_CPU_code();

Pour de plus amples informations voir ce pdf.

3.2.2. Exemple

On considère l'exemple de la figure suivante pour laquelle on dispose :

• d'une grille 2D de $(3 × 4)$ blocs,
• chaque bloc 2D est composé de $(2 × 5)$ threads

CUDA définit pour le kernel à exécuter des variables de type dim3 :

• gridDim : dimensions de la grille (3, 4, 1)
• blockDim : dimensions d'un bloc (2, 5, 1)
• blockIdx : coordonnées du bloc courant à l'intérieur de la grille
• threadIdx : coordonnées du thread courant à l'intérieur d'un bloc

Au final on dispose d'un ensemble de $(3 × 4) × (2 × 5) = 120$ threads à exécuter qui seront exécutés par blocs de $(2 × 5)$ threads. A l'intérieur du kernel, pour connaître les coordonnées globales d'un thread par rapport à la grille, on utilise les formules suivantes (si la grille est de type 2D et si un bloc est de type 2D) :

int x = blockIdx.x * blockDim.x + threadIdx.x
int y = blockIdx.y * blockDim.y + threadIdx.y

Si chaque thread correspond à un élément d'un tableau stocké de manière linéaire en mémoire, on détermine l'offset de l'élément par la formule :

int offset = y * gridDim.x * blockDim.x + x

Par exemple, le thread (y=1,x=3) du bloc (y=1,x=2) a pour valeurs :

x = 2 * 5 + 3 = 13
y = 1 * 2 + 1 = 3
offset = 3 * 4 * 5 + 13 = 73

Comme on peut le voir sur l'exemple précédent, les threads dans un bloc accèdent des éléments d'un tableau de manière linéaire pour chaque ligne. On verra qu'il est plus intéressant d'accèder des éléments consécutifs pour chaque bloc de threads, sinon on risque de dégrader les performances.

3.3. Etude de cas

3.3.1. Somme de longs vecteurs

Considérons le cas de la somme de longs vecteurs 1D. Quel doit être le choix de la taille de la grille ou la taille d'un bloc ?

3.3.1.a  solution fixe avec un seul bloc + stride

Un première solution consiste à utiliser une grille composée d'un seul bloc, le bloc étant composé de plusieurs threads et on itère jusqu'à atteindre la taille du vecteur en se déplaçant du nombre de threads (=stride, en français foulée, enjambée, avancée).

const int VECTOR_SIZE = 20480; // = 512 * 40
const int NUMBER_OF_THREADS = 512;

__global__ void kernel_sum(float *x, float *y, float *z, int size) {
	int tid = threadIdx.x;
	while (tid < size) {
		z[tid] = x[tid] + y[tid];
		// stride of blockDim.x
		tid += blockDim.x;
	}
}

kernel_sum<<< 1, NUMBER_OF_THREADS >>>(x_d, y_d, z_d, VECTOR_SIZE); 

Pour cette solution qui n'utilise que 512 threads:

• le thread 0 traitera les éléments du tableau d'indice : 0, 512, 1024, ...
• le thread 1 traitera les éléments du tableau d'indice : 1, 513, 1025, ...
• le thread i traitera les éléments du tableau d'indice : i, 512+i, 1024+i, ...
• le thread 511 traitera les éléments du tableau d'indice : 511, 1023, 1535, ...

On note que pour cette solution l'occupancy est faible, il en résulte des temps de calcul élevés par rapport à la solution suivante.

3.3.1.b  solution adaptative avec plusieurs blocs

Une deuxième approche consiste à adapter la taille de la grille ou du bloc à la taille du vecteur en utilisant autant de blocs que nécessaire. Dans ce cas le nombre total de threads donnés par la taille de la grille et du bloc doit être supérieur ou égal au nombre de données à traiter.

Pour un vecteur de taille 20480 (20 * 1024) floats :

• si la taille maximale d'un bloc (sur la dimension x) est de 512, alors on a 40 (= 20480 / 512) blocs dans la grille
• si la taille maximale d'un bloc (sur la dimension x) est de 1024, alors on a 20 blocs dans la grille

const int VECTOR_SIZE = 20480;
const int NUMBER_OF_THREADS = 512;
int NUMBER_OF_BLOCKS = (VECTOR_SIZE + NUMBER_OF_THREADS - 1) / NUMBER_OF_THREADS;

__global__ void kernel_sum(float *x, float *y, float *z, int size) {
	int gtid = threadIdx.x + blockIdx.x * blockDim.x;
	z[gtid] = x[gtid] + y[gtid];
}

kernel_sum<<< NUMBER_OF_BLOCKS, NUMBER_OF_THREADS >>>(x_d, y_d, z_d, VECTOR_SIZE); 

Quelle solution est la plus intéressante ? Vaut-il mieux avoir peu ou beaucoup de blocs ? Dans le cas présent, on pourrait choisir de prendre des blocs de 32 threads, ce qui donnerait une grille de 640 blocs.

Que faire si la taille de la grille n'est pas un multiple de 16, 32, ... ?

Voici quelques résultats sur GeForce GTX 560 Ti (10.000 itérations pour vecteur de taille 262.144 floats) :

On note que les temps les plus intéressants correspondent à une occupancy généralement maximum mais pas forcément.

3.3.1.c  Nombre de threads important mais inférieurs au nombre de données à traiter

Enfin une troisième solution, similaire à la première consiste à disposer de nombreux threads (plus de 1024) mais moins que de données à traiter et à réaliser un stride.

const int VECTOR_SIZE = 20480;
const int NUMBER_OF_THREADS = 512;
int NUMBER_OF_BLOCKS = 4;

__global__ void kernel_sum(float *x, float *y, float *z, int size) {
	int gtid = threadIdx.x + blockIdx.x * blockDim.x;
	int stride = gridDim.x * blockDim.x;
	 	
	while (gtid < size) {
	   z[gtid] = x[gtid] + y[gtid];
	   gtid += stride;
	}   
}

kernel_sum<<< NUMBER_OF_BLOCKS, NUMBER_OF_THREADS >>>(x_d, y_d, z_d, VECTOR_SIZE); 

Dans cet exemple on génère $4 × 512 = 2048$ threads, le stride est donc de 2048, on aura donc 10 itérations de la boucle while.

3.3.2. Mesure des performances

3.3.2.a  Mesures de temps

Les mesures de performance des deux exemples précédents sont réalisées grâce à l'API CUDA.

CUDA fournit des compteurs ainsi que des fonctions associées pour gérer ces compteurs.

L'appel à cudaEventCreate réalise une allocation mémoire d'une structure de données, il faudra donc libérer l'espace alloué grâce à cudaEventDestroy.

3.3.2.b  Profilage : anciennes versions avec nvprof

On peut également utiliser l'utilitaire nvprof :

nvprof --unified-memory-profiling off program arguments

ou son équivalent graphique nvvp (NVidia Visual Profiler).

Penser à utiliser également --print-gpu-trace ou --metrics achieved_occupancy

3.3.2.c  Profilage : à partir de Compute Capability 8.0

nvprof n'est plus supporté sur les GPU possédant une Compute capability 8.0 et plus.

Il est conseillé d'utiliser NVidia NSight ou nsys.

Utiliser les commandes suivantes :

nsys profile -o profile_report --force-overwrite true build/bin/vec_sum.exe 1000000
WARNING: CPU IP/backtrace sampling not supported, disabling.
Try the 'nsys status --environment' command to learn more.

WARNING: CPU context switch tracing not supported, disabling.
Try the 'nsys status --environment' command to learn more.

Collecting data...
total.cpu=30000
total.gpu=30000
Generating '/tmp/nsys-report-b083.qdstrm'
[1/1] [========================100%] profile_report.nsys-rep
Generated:
    .../profile_report.nsys-rep

Cette commande (nsys profile) génère un fichier profile_report.nsys-rep qui peut ensuite être analysé :

• cuda_api_summary : indique les fonctions autres que les kernels utilisées (API)
• cuda_gpu_mem_time_summary : donne des informations sur le temps mis pour réaliser les transferts mémoire
• cuda_gpu_kern_summary : indique les kernels utilisées
• cuda_gpu_mem_size_summary : donne des informations sur le nombre d'octets liés aux transferts mémoire

\$ nsys stats -r cuda_api_sum,cuda_gpu_mem_time_sum,cuda_gpu_kern_sum,cuda_gpu_mem_size_sum --force-export=true profile_report.nsys-rep
Generating SQLite file profile_report.sqlite from profile_report.nsys-rep
Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/reports/cuda_api_sum.py]... 

 ** CUDA API Summary (cuda_api_sum):

 Time (%)  Total Time (ns)  Num Calls    Avg (ns)     Med (ns)    Min (ns)    Max (ns)   StdDev (ns)            Name         
 --------  ---------------  ---------  ------------  -----------  ---------  ----------  ------------  ----------------------
     92,0       60 506 544          3  20 168 848,0     72 586,0     64 751  60 369 207  34 814 532,0  cudaMalloc            
      5,0        3 776 070          3   1 258 690,0  1 290 426,0  1 105 370   1 380 274     140 172,0  cudaMemcpy            
      1,0          647 282          3     215 760,0    258 534,0    111 118     277 630      91 124,0  cudaFree              
      0,0          285 325          1     285 325,0    285 325,0    285 325     285 325           0,0  cudaLaunchKernel      
      0,0              521          1         521,0        521,0        521         521           0,0  cuModuleGetLoadingMode

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/reports/cuda_gpu_mem_time_sum.py]... 

 ** CUDA GPU MemOps Summary (by Time) (cuda_gpu_mem_time_sum):

 Time (%)  Total Time (ns)  Count   Avg (ns)     Med (ns)    Min (ns)   Max (ns)   StdDev (ns)           Operation          
 --------  ---------------  -----  -----------  -----------  ---------  ---------  -----------  ----------------------------
     67,0        2 426 704      2  1 213 352,0  1 213 352,0  1 213 256  1 213 448        135,0  [CUDA memcpy Host-to-Device]
     33,0        1 193 255      1  1 193 255,0  1 193 255,0  1 193 255  1 193 255          0,0  [CUDA memcpy Device-to-Host]

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/reports/cuda_gpu_kern_sum.py]... 

 ** CUDA GPU Kernel Summary (cuda_gpu_kern_sum):

 Time (%)  Total Time (ns)  Instances  Avg (ns)   Med (ns)   Min (ns)  Max (ns)  StdDev (ns)                      Name                     
 --------  ---------------  ---------  ---------  ---------  --------  --------  -----------  ---------------------------------------------
    100,0          127 745          1  127 745,0  127 745,0   127 745   127 745          0,0  kernel_sum(double *, double *, double *, int)

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/reports/cuda_gpu_mem_size_sum.py]... 

 ** CUDA GPU MemOps Summary (by Size) (cuda_gpu_mem_size_sum):

 Total (MB)  Count  Avg (MB)  Med (MB)  Min (MB)  Max (MB)  StdDev (MB)           Operation          
 ----------  -----  --------  --------  --------  --------  -----------  ----------------------------
 16,000          2  8,000     8,000     8,000     8,000     0,000        [CUDA memcpy Host-to-Device]
 8,000           1  8,000     8,000     8,000     8,000     0,000        [CUDA memcpy Device-to-Host]

On voit que :

• cudaMalloc a été appelée 3 fois pour allouer les vecteurs gpu_vec_x, gpu_vec_y et gpu_vec_z
• le temps total pour l'allocation des trois vecteurs est de 60 ms
• la copie des données de cpu_vec_x et cpu_vec_y vers gpu_vec_x et gpu_vec_y est de 2 ms
• la copie des données de gpu_vec_z vers cpu_vec_z a duré 1,19 ms
• l'exécution du kernel a pris 127 micro seconds
• un vecteur d'un million de doubles, soit $8 × 10^6$, occupe $8 × 10^6 / (1024^2) = 7,62$ Mo mais CUDA reporte 8 Mo

Note : la mémoire allouée est généralement multiple de blocks de 4 ko de manière à avoir des données coalescentes (soudées).

Note : pour obtenir les temps en milli secondes il faut utiliser --timeunit milli.

\$ nsys analyze --force-export=true profile_report.nsys-rep
Generating SQLite file profile_report.sqlite from profile_report.nsys-rep
Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/cuda_memcpy_async.py]... 

 ** CUDA Async Memcpy with Pageable Memory (cuda_memcpy_async):

There were no problems detected related to memcpy operations using pageable
memory.

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/cuda_memcpy_sync.py]... 

 ** CUDA Synchronous Memcpy (cuda_memcpy_sync):

The following are synchronous memory transfers that block the host. This does
not include host to device transfers of a memory block of 64 KB or less.

Suggestion: Use cudaMemcpy*Async() APIs instead.

 Duration (ns)  Start (ns)   Src Kind  Dst Kind  Bytes (MB)   PID    Device ID  Context ID  Green Context ID  Stream ID      API Name    
 -------------  -----------  --------  --------  ----------  ------  ---------  ----------  ----------------  ---------  ----------------
     1 213 448  333 410 666  Pageable  Device    8,000       36 549          0           1                            7  cudaMemcpy_v3020
     1 213 256  334 703 314  Pageable  Device    8,000       36 549          0           1                            7  cudaMemcpy_v3020
     1 193 255  336 132 955  Device    Pageable  8,000       36 549          0           1                            7  cudaMemcpy_v3020

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/cuda_memset_sync.py]... 
SKIPPED: profile_report.sqlite could not be analyzed because it does not contain the required CUDA data. Does the application use CUDA memset APIs?

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/cuda_api_sync.py]... 

 ** CUDA Synchronization APIs (cuda_api_sync):

There were no problems detected related to synchronization APIs.

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/gpu_gaps.py]... 

 ** GPU Gaps (gpu_gaps):

There were no problems detected with GPU utilization. GPU was not found to be
idle for more than 500ms.

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/gpu_time_util.py]... 

 ** GPU Time Utilization (gpu_time_util):

The following are time regions with an average GPU utilization below 50%%.
Addressing the gaps might improve application performance.

Suggestions:
   1. Use CPU sampling data, OS Runtime blocked state backtraces, and/or OS
Runtime APIs related to thread synchronization to understand if a sluggish or
blocked CPU is causing the gaps.
   2. Add NVTX annotations to CPU code to understand the reason behind the gaps.

 Row#  In-Use (%)  Duration (ns)  Start (ns)    PID    Device ID  Context ID
 ----  ----------  -------------  -----------  ------  ---------  ----------
    1        37,0        130 519  335 890 527  36 549          0           1
    2        39,0        130 519  334 585 337  36 549          0           1

Processing [profile_report.sqlite] with [/opt/nvidia/nsight-systems/2024.5.1/host-linux-x64/rules/dx12_mem_ops.py]... 
SKIPPED: profile_report.sqlite could not be analyzed because it does not contain the required DX12 data. Does the application use DX12 APIs?