Ecole Doctorale MaSTIC

26 et 27 Juin 2024, Salle H113
Faculté des Sciences
Université d'Angers

1.1. Cours

Vous trouverez une partie du cours ici

1.2. Ressources

1.3. Mise en application

Les travaux pratiques ont lieu en salle H113 qui comprend des machines dotées des caractéristiques suivantes :

1.3.1. Somme de vecteurs

Il s'agit d'un exemple de base, classique, qui permet de comprendre l'enchaînement des opérations et l'interaction entre CPU et GPU.

On dispose de trois vecteurs $x$, $y$ et $z$ de taille $N$ et on calcule :

$∀i ∈ [0..N[, z[i] = x[i] + y[i]$.

1.3.1.a  Programme

Créer un programme appelé vec_sum.cu qui réalisera le calcul. Ce programme fera la somme de vecteurs de float ou de double et prendra en paramètre la taille des vecteurs. On pourra par exemple l'exécuter avec la commande :

./vec_sum.exe -n 1000000

Vous pouvez récupérer l'archive cuda_ed_etud.tgz qui contient quelques fichiers nécessaires à l'élaboration des programmes que nous allons mettre en oeuvre.

1.3.1.b  Vecteurs sur le CPU

On définira le type real (comme float ou double) puis on déclare les vecteurs sur le CPU :

  1. #include <cstdint>
  2. #include <iostream>
  3. #include <iomanip>
  4. using namespace std;
  5.  
  6. typedef float real;
  7.  
  8. // équivalent du type size_t
  9. typedef uint32_t u32;
  10.  
  11. // size of a vector as the number of elements
  12. u32 vec_size;
  13.  
  14. // pointers to vectors that will be allocated
  15. real *cpu_vec_x;
  16. real *cpu_vec_y;
  17. real *cpu_vec_z;

Un fois qu'on a déterminé la taille des vecteurs (paramètre du programme), on peut allouer ceux-ci :

  1. cpu_vec_x = new real[vec_size];
  2. cpu_vec_y = new real[vec_size];
  3. cpu_vec_z = new real[vec_size];

puis les instancier :

  1. for (u32 i = 0; i < vec_size; ++i)
  2. {
  3.     cpu_vec_x[i] = 0.1;
  4.     cpu_vec_y[i] = 0.2;
  5. }

1.3.1.c  Vecteurs sur le GPU

On commence par déclarer les vecteurs puis on les alloue sur le GPU :

  1. real *gpu_vec_x;
  2. real *gpu_vec_y;
  3. real *gpu_vec_z;
  4.  
  5.  
  6. cudaMalloc( (void **) &gpu_vec_x, vec_size * sizeof(real) );
  7. cudaMalloc( (void **) &gpu_vec_y, vec_size * sizeof(real) );
  8. cudaMalloc( (void **) &gpu_vec_z, vec_size * sizeof(real) );

Il faut ensuite copier les données des vecteurs $x$ et $y$ vers les vecteurs alloués sur le GPU :

  1. // cudaMemcpy(destination, source, taille en octets, sens de la copie)
  2. cudaMemcpy( gpu_vec_x, cpu_vec_x, vec_size* sizeof(real), cudaMemcpyHostToDevice );
  3. cudaMemcpy( gpu_vec_y, cpu_vec_y, vec_size* sizeof(real), cudaMemcpyHostToDevice );

Avant de lancer le kernel qui réalisera le calcul, on calcule le nombre de blocs en fonction de la taille d'un bloc. Dans le cas présent on prend un bloc de taille maximale c'est à dire de 1024 threads.

  1. int block_dim = 1024;
  2. int nbr_blocks = (vec_size + block_dim - 1) / block_dim;

On lance le kernel de calcul sur le GPU :

  1. kernel_sum<<<nbr_blocks, block_dim>>>( gpu_vec_x, gpu_vec_y, gpu_vec_z, vec_size );

Enfin on copie le résultat qui se trouve dans gpu_vec_z vers cpu_vec_z :

  1. // cudaMemcpy(destination, source, taille en octets, sens de la copie)
  2. cudaMemcpy( cpu_vec_z, gpu_vec_z, vec_size* sizeof(real), cudaMemcpyDeviceToHost );

On affiche le vecteur obtenu :

  1. for (u32 i = 0; i < vec_size; ++i) {
  2.     cout << i << " " << cpu_vec_z[i] << endl;
  3. }

On libère la mémoire :

  1. delete [] cpu_vec_x;
  2. delete [] cpu_vec_y;
  3. delete [] cpu_vec_z;
  4.  
  5. cudaFree( gpu_vec_x );
  6. cudaFree( gpu_vec_y );
  7. cudaFree( gpu_vec_z );
  8.  

1.3.1.d  Le kernel

On utilise autant de threads que d'éléments dans le vecteur, le kernel s'écrit donc :

  1. __global__
  2. void kernel_sum( real *x, real *y, real *z, int size)
  3. {
  4.     int tid = blockDim.x * blockIdx.x + threadIdx.x;
  5.  
  6.     if (tid < size) {
  7.         z[tid] = x[tid] + y[tid];
  8.     }
  9.  
  10. }

1.3.1.e  Compilation

On procède en deux étapes en utilisant nvcc le compilateur CUDA :

  1. nvcc --compile -o vec_sum.o vec_sum.cu  
  2.     -gencode arch=compute_52,code=sm_52
  3.     -gencode arch=compute_60,code=sm_60  
  4.     -gencode arch=compute_61,code=sm_61  
  5.     -gencode arch=compute_75,code=sm_75  
  6.     -gencode arch=compute_80,code=sm_80  
  7.     -gencode arch=compute_86,code=sm_86  
  8.     -std=c++11 -g -G -Xcompiler="-std=c++11 -O3 -Wall"
  9.  
  10.  
  11. nvcc --link -o vec_sum.exe vec_sum.o  

Dans cet exemple on a spécifié plusieurs Compute Capability (CC = 5.2, 6.0, 6.1, 7.5, 8.0, 8.6) car on ne connaît pas le type de GPU utilisé. En général, on ne mettra qu'une seule CC, celle du GPU que l'on utilisera.

1.3.1.f  Tests

Réaliser quelques tests avec des tailles de vecteurs de 1.000 puis 1.000.000 d'éléments.

Que remarque t-on ?

1.3.2. Réduction

1.3.2.a  Principe

La réduction consiste à calculer, par exemple, la somme des éléments d'un vecteur :

$∀i ∈ [0..N[,\; sum = ∑↙{i=0}↖{n-1} z[i]$

La réduction est un algorithme problématique mais pour lequel on trouve rapidement des solutions qui permettent de conserver le parallèlisme.

1.3.2.b  Sur le CPU

Sur un CPU doté de $k$ threads de calcul, on divisera le calcul de la réduction en $k$ sommes partielles, puis un seul thread réalisera la somme des sommes partielles.

  1. float sum = 0.0;
  2.  
  3. #pragma omp for reduction(+ : sum)
  4. for (int i = 0; i < vec_size; ++i)
  5. {
  6.     sum += z[i];
  7. }

1.3.2.c  Première implantation sur le GPU

On peut reprendre ce principe sur GPU en utilisant un seul bloc de threads :

  1. __global__
  2. void kernel_reduction_1(real *z, int vec_size)
  3. {
  4.     int tid = threadIdx.x;
  5.     int i = tid + blockDim.x;
  6.  
  7.     while (i < vec_size) {
  8.         z[tid] += z[i];
  9.         i += blockDim.x;
  10.     }
  11.  
  12. }

Par exemple si on utilise un bloc de 1024 threads et des vecteurs de taille 10240, on parcourera 9 blocs, le premier bloc contiendra les sommes partielles.

On transfère le premier bloc vers le CPU qui réalise la somme des sommes partielles.

1.3.2.d  Deuxième implantation sur le GPU

On va utiliser de la mémoire partagée, celle-ci agit comme une mémoire cache et son temps accès est de l'ordre de 10 à 20 cycles d'horloge, tandis que la mémoire GDDR peut avoir un temps d'accès de l'ordre de centaines de cycles d'horloge.

  1. __global__
  2. void kernel_reduction_2(real *z, int size) {
  3.    
  4.     // déclaration de la mémoire partagée
  5.     // la taille n'est pas précisée, elle le sera
  6.     // lors de l'appel du kernel
  7.     extern __shared__ real shared_mem[];
  8.  
  9.     // thread index
  10.     int tid = threadIdx.x;
  11.     shared_mem[tid] = z[tid];
  12.  
  13.     // index;
  14.     int i = tid + blockDim.x;
  15.  
  16.     while (i < size)
  17.     {
  18.         shared_mem[tid] += z[i];
  19.         i += blockDim.x;
  20.     }
  21.  
  22.     z[tid] = shared_mem[tid];
  23. }
  24.  
  25. // l'appel au kernel se fait en précisant la taille de la mémoire partagée
  26. // en troisième paramètre ici un tableau de 1024 real, si block_size=1024
  27. kernel_reduc<<<1, block_size, block_size*sizeof(real)>>>(gpu_vec_z, vec_size);

L'influence de la mémoire partagée sera limitée puisqu'on accède principalement à la mémoire du GPU (GDDR).

Il faudra comme précédemment transférer le premier bloc vers la mémoire du CPU pour réaliser la somme des sommes partielles.

1.3.2.e  Troisième implantation sur le GPU

Dans cette dernière approche qui est sensée être la plus rapide, la somme des sommes partielles est réalisée sur le GPU et est placée dans $z[0]$ :

  1. __global__
  2. void kernel_reduc(real *z, int size) {
  3.    
  4.     extern __shared__ real shared_mem[];
  5.  
  6.     // thread index
  7.     int tid = threadIdx.x;
  8.     shared_mem[tid] = z[tid];
  9.  
  10.     // index;
  11.     int i = tid + blockDim.x;
  12.  
  13.     while (i < size)
  14.     {
  15.         shared_mem[tid] += z[i];
  16.         i += blockDim.x;
  17.     }
  18.  
  19.     // calcul de la somme des sommes partielles sur le GPU
  20.     __syncthreads();
  21.    
  22.     int stride = blockDim.x / 2;
  23.    
  24.     while (stride > 0) {
  25.         if (tid < stride) {
  26.             shared_mem[tid] += shared_mem[tid + stride];
  27.         }
  28.         stride = stride / 2;
  29.         __syncthreads();
  30.     }
  31.  
  32.     if (tid == 0) {
  33.         z[0] = shared_mem[0];
  34.     }
  35.  
  36. }
  37.  

Ici il est nécessaire de synchroniser les calculs à chaque étape :

réduction sur GPU

A chaque étape on divise le bloc en deux et on additionne la partie haute du bloc à la partie basse, puis on recommence.

La fonction __syncthreads() permet :

  • la synchronisation de tous les threads dans un bloc donné. Lorsque cette fonction est appelée, tous les threads doivent atteindre cette instruction avant que l'exécution ne puisse continuer. Cela signifie que tous les threads du bloc doivent terminer leurs calculs jusqu'à ce point avant que tout thread puisse avancer.
  • d'assurer la consistance de la mémoire partagée : toutes les écritures dans la mémoire partagée faites par les threads sont visibles par tous les autres threads du bloc. Cela est crucial lorsque les threads écrivent dans la mémoire partagée, puis lisent ces valeurs. Sans synchronisation, certains threads pourraient lire des valeurs non mises à jour.

1.4. Courbes de Julia

Voir ce TP.

Il s'agit d'un exemple typique de parallélisme ou le GPU peut se révéler plus intéressant que le CPU.

1.5. CUDA avec Python

Installer le module pycuda.

Voici un exemple avec la somme de vecteurs :

  1. import pycuda.autoinit
  2. import pycuda.driver as cuda
  3. import numpy as np
  4. from pycuda.compiler import SourceModule
  5.  
  6. # Kernel CUDA pour additionner les vecteurs
  7. kernel_code = """
  8. __global__ void vector_sum(float *x, float *y, float *z, int size)
  9. {
  10.    int tid = threadIdx.x + blockIdx.x * blockDim.x;
  11.    if (tid < size)
  12.    {
  13.        z[tid] = x[tid] + y[tid];
  14.    }
  15. }
  16. """
  17.  
  18. # Taille des vecteurs
  19. n = 1024
  20. # Créer des vecteurs aléatoires
  21. a = np.random.randn(n).astype(np.float32)
  22. b = np.random.randn(n).astype(np.float32)
  23. c = np.zeros_like(a)
  24.  
  25. # Allouer de la mémoire sur le GPU et copier les données
  26. a_gpu = cuda.mem_alloc(a.nbytes)
  27. b_gpu = cuda.mem_alloc(b.nbytes)
  28. c_gpu = cuda.mem_alloc(c.nbytes)
  29.  
  30. cuda.memcpy_htod(a_gpu, a)
  31. cuda.memcpy_htod(b_gpu, b)
  32.  
  33. # Compiler le kernel
  34. mod = SourceModule(kernel_code)
  35. vector_sum_kernel = mod.get_function("vector_sum")
  36.  
  37. # Définir la taille du bloc et du grid
  38. block_size = 256
  39. grid_size = (n + block_size - 1) // block_size
  40.  
  41. # Lancer le kernel
  42. vector_sum_kernel(a_gpu, b_gpu, c_gpu, np.int32(n), block=(block_size, 1, 1), grid=(grid_size, 1, 1))
  43.  
  44. # Copier le résultat du GPU vers le CPU
  45. cuda.memcpy_dtoh(c, c_gpu)
  46.  
  47. # Vérifier le résultat
  48. print("Vecteur a:", a)
  49. print("Vecteur b:", b)
  50. print("Somme a + b:", c)
  51.  
  52.  

1.6. Analyse

Comparer les caractéristiques des cartes suivantes :

On s'intéressera notamment aux caractéristiques suivantes :

Combien de RTX 4090 faut-il pour égaler la puissance d'une A100 ?

1.7. Vidéo

Vidéo