CUDA : cours    travaux pratiques

Sommaire

6. TD Réduction

6.1. Introduction

L'algorithme de réduction est un algorithme important en parallèlisme.

6.1.1. Réduction - définition

On considère un tableau de $n$ éléments et un opérateur associatif et commutatif $⊕$ doté d'un élément neutre $I$. A partir du tableau

$$ [ a_0, a_1, ..., a_{n-1} ] $$

on désire obtenir un unique résultat :

$$ a_0 ⊕ a_1 ⊕ ... ⊕ a_{n-1} $$

Si $⊕$ est l'addition alors la réduction permet d'obtenir la somme des éléments du tableau.

On peut appliquer cet algorithme avec d'autres opérateurs comme la multiplication, minimum, maximum.

Pour que l'algorithme soit parallélisable efficacement il suffit que l'opérateur $⊕$ soit associatif et commutatif :

6.2. Version séquentielle

Une version séquentielle est donnée par le code EZ language suivant : on calcule la somme.

function reduction(input is array of integer) return integer
local i, sum are integer 
begin
	sum = 0
	for i in input.range() do
		sum = sum + input[i];
	endfor
	return sum
end reduction

Voici une version plus élaborée : on utilise la STL et les template afin de simplifier l'écriture du code.

On utilise un tableau avec les données en entrée input.

  1. // compile with g++ -std=c++11 -Wall -O2
  2. #include <iostream>
  3. #include <iterator>
  4. #include <algorithm>
  5. #include <cstdlib>
  6. #include <cmath>
  7. #include <cstring>
  8. using namespace std;
  9.  
  10. /**
  11.  * print array
  12.  */
  13. template<class T>
  14. void print_array(string s, T *a, int size) {
  15.     cout << s;
  16.     copy(&a[0], &a[size-1], ostream_iterator<T>(cout, ", "));
  17.     cout << a[size-1];
  18.     cout << endl;
  19. }
  20.  
  21. /**
  22.  * functor for addition
  23.  */
  24. template<class T>
  25. class PlusOperator {
  26. public:
  27.     static T operation(T a, T b) {
  28.         return a + b;
  29.     }
  30.    
  31.     static T identity() {
  32.         return 0;
  33.     }
  34. };
  35.  
  36. /**
  37.  * functor for product
  38.  */
  39. template<class T>
  40. class ProductOperator {
  41. public:
  42.     static T operation(T a, T b) {
  43.         return a * b;
  44.     }
  45.    
  46.     static T identity() {
  47.         return 1;
  48.     }
  49. };
  50.  
  51. /**
  52.  * Reduction with functor as template parameter
  53.  * @param in  input array
  54.  * @param size number of elements of the array
  55.  */
  56. template<class T, class Oper>
  57. T reduction(T *input, int size) {
  58.     // set first output element as first input element
  59.     T result = Oper::identity();
  60.     for (int i=0; i<size; ++i) {
  61.         result = Oper::operation(result, input[i]);
  62.     }
  63.     return result;
  64. }
  65.  
  66.  
  67. /**
  68.  * main function
  69.  */
  70. int main(int argc, char *argv[]) {
  71.     int MAX_ELEMENTS = 16;
  72.     int *int_input;
  73.     double *dbl_input;
  74.     int start = 1;
  75.     int int_result;
  76.     double dbl_result;
  77.    
  78.     if (argc > 1) {
  79.         MAX_ELEMENTS = atoi(argv[1]);
  80.     }
  81.    
  82.     // allocate input and output arrays
  83.     int_input = new int [ MAX_ELEMENTS ];
  84.     dbl_input = new double [ MAX_ELEMENTS ];
  85.    
  86.     // fill int_input array with 1, 2, 3, ....
  87.     generate(&int_input[0], &int_input[MAX_ELEMENTS], [&start]() {
  88.         return start++;
  89.     });
  90.  
  91.     // fill dbl_input
  92.     start = 1;
  93.     generate(&dbl_input[0], &dbl_input[MAX_ELEMENTS], [&start]() {
  94.         return start++;
  95.     });
  96.  
  97.     // print input values
  98.     print_array<int>("input array of integers : ", int_input, MAX_ELEMENTS);
  99.     print_array<double>("input array of integers : ", dbl_input, MAX_ELEMENTS);
  100.    
  101.     // reduction with integer and +
  102.     int_result = reduction<int, PlusOperator<int>>(int_input, MAX_ELEMENTS);
  103.     cout << "result with + : " << int_result << endl;
  104.     // reduction with double and *
  105.     dbl_result = reduction<double, ProductOperator<double>>(dbl_input, MAX_ELEMENTS);
  106.     cout << "result with * : " << fixed << dbl_result << endl;
  107.    
  108.     delete [] int_input;
  109.     delete [] dbl_input;
  110.    
  111.     return EXIT_SUCCESS;
  112. }
  113.  
  114.  
./code_reduction_seq.exe 10
input array of integers : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
input array of integers : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
result with + : 55
result with * : 3628800.000000

La complexité de cet algorithme est en $O(n)$, ou $n$ est la taille du tableau.

6.3. Versions parallèles

6.3.1. le coeur de l'algorithme

On peut écrire une version qui fonctionne comme la version séquentielle mais en utilisant une addition atomique :

function parallel_reduction(input is array of integer) return integer
local i, sum are integer 
begin
	sum = 0
	meta parallel
	for i in input.range() do
		sum = atomic_add(sum, input[i]);
	endfor 
	return sum
end reduction

La fonction atomic_add est bloquante et seul un thread peut l'exécuter à la fois. On obtient au final un traitement séquentiel, ce qui n'est pas intéressant du point de vue de la performance et nous fait perdre l'utilité du parallèlisme.

L'algorithme dans sa version parallèle efficace utilise une décomposition arborescente du tableau.

Voici un exemple avec l'addition :

Reduction

Ainsi la complexité est en $O(\log_2(n))$ : on effectue toujours $n-1$ additions mais en $\log_2(n)$ étapes. Dans le cas présent :

6.3.1.a  version 0

L'algorithme peut être écrit sous la forme suivante. On considère un tableau input composé de $n$ éléments, le premier élement du tableau étant input[0] et le dernier input[n-1].

On considère ici que $n$ est une puissance de 2 supérieure ou égale à 1024.

On utilise une grille 1D de blocs de 256 threads par exemple de manière à avoir au moins 4 blocs (puisque $n ≧ 1024)$.

On parcourt le tableau de deux en deux, puis de quatre en quatre, ..., grâce à la variable stride que l'on peut traduire en français par foulée ou enjambée.

  1. /**
  2.  * kernel to perform reduction of array used as
  3.  * input and output
  4.  */
  5. __global__ void reduction_v0(T *input, int size) {
  6.     // we use a grid 1D of blocks 1D
  7.     // Global thread Index
  8.     int i = blockIdx.x * blockDim.x + threadIdx.x;
  9.  
  10.     // perform reduction in GPU global memory
  11.     for (int stride = 1; stride < size; stride *= 2) {
  12.         if (((i + 1) % (2*stride)) == 0) {
  13.             input[i] += input[i - stride];
  14.         }
  15.         __syncthreads();
  16.     }
  17. }
  18.  
  19. // call to kernel
  20. // consider that size = 128
  21. dim3 block_dim(64, 1, 1);
  22. dim3 grid_dim(2, 1, 1);
  23. reduction_v0<<<grid_dim, block_dim>>>(gpu_input, size);
  24.  
  25.  
1, 2, 3,  4, 5,  6, 7,  8  
1, 3, 3,  7, 5, 11, 7, 15
1, 3, 3, 10, 5, 11, 7, 26
1, 3, 3, 10, 5, 11, 7, 36

On peut modifier le code pour obtenir une autre version qui commence avec les indices pairs:

  1. /**
  2.  * kernel to perform reduction of array used as
  3.  * input and output
  4.  */
  5. __global__ void reduction_v0(T *input, int size) {
  6.     // we use a grid 1D of blocks 1D
  7.     // Global thread Index
  8.     int i = blockIdx.x * blockDim.x + threadIdx.x;
  9.  
  10.     // perform reduction in GPU global memory
  11.     for (int stride = 1; stride < size; stride *= 2) {
  12.         if ((i % (2*stride)) == 0) {
  13.             input[i] += input[i + stride];
  14.         }
  15.         __syncthreads();
  16.     }
  17. }
  18.  
  19. // call to kernel
  20. // consider that size = 128
  21. dim3 block_dim(64, 1, 1);
  22. dim3 grid_dim(2, 1, 1);
  23. reduction_v0<<<grid_dim, block_dim>>>(gpu_input, size);
  24.  
  25.  
 1, 2, 3, 4,  5, 6,  7, 8
 3, 2, 7, 4, 11, 6, 15, 8
10, 2, 7, 4, 26, 6, 15, 8
36, 2, 7, 4, 26, 6, 15, 8

Analyse

  • on modifie les données en entrée : on pourrait vouloir les garder intactes
  • on travaille en mémoire globale du GPU : input[i] += input[i - stride]; ce qui n'est pas très performant
  • problème majeur : cela ne fonctionne pas ! En effet sous GTX 560 Ti avec un tableau initial composé de '1' et l'opérateur +, l'algorithme dans sa version v0_2 semble fonctionner jusqu'à size=8192, mais donne 8193 pour size=16384 alors qu'on devrait obtenir 16384.

    Cela est du au fait qu'il n'y a pas de synchronisation des écritures en mémoire globale du GPU en fonction des blocs.

    Si on prend l'algorithme dans sa version v0_1 alors cela fonctionne jusqu'à size=32768 mais à partir de size=65536 on obtient parfois le résultat escompté et parfois un résultat différent (61440 par exemple)
  • l'instruction __syncthreads() permet uniquement de synchroniser les threads d'un bloc, mais il n'existe pas d'instruction de synchronisation de tous les blocs ou entre les blocs.
  • l'explication du non fonctionnement est la suivante (Forum NVidia) :
    I am to use elementary numbers, to make it easier; the principle applies
    nevertheless.
    You are using device x, with y sm's, with the constraints and conditions of
    device x's compute capability you wish to scan/reduce 2000 elements your device
    can only seat 10 blocks of 100 threads each, at a time, meaning 1000 elements

    input[i] += input[i + stride]

    when thread i gets to a stride that would imply i + stride = 1500, it expects
    element 1500 is ready but element 1500 can not be ready; in fact, element 1001 to
    2000 would not be ready and so the reduction/scan breaks down, because of a violation
    of its premises

Améliorations

  • calculer une somme partielle pour chaque bloc et stocker les résultats dans un autre tableau
  • charger les données en mémoire partagée et travailler en mémoire partagée : accélération des traitements et synchronisation des données/threads

6.3.1.b  version 1

Le principe de la version 1 est le suivant :

Réduction version 1

  1. /**
  2.  * reduction
  3.  * @param input array of input data
  4.  * @param output partial sum of each block
  5.  * @param size number of elements of ind
  6.  */
  7. __global__ void reduction_v1(int *input, int *output, int size) {
  8.     // local thread index
  9.     int ltid = threadIdx.x;
  10.     // global thread index
  11.     int gtid = blockIdx.x * blockDim.x + threadIdx.x;
  12.  
  13.     extern __shared__ int smem[];
  14.    
  15.     // load data in shared memory
  16.     smem[ltid] = (i < size) ? input[gtid] : 0;
  17.     __syncthreads();
  18.  
  19.     // perform reduction in shared memory
  20.     for (int stride = 1; stride < blockDim.x; stride *= 2) {
  21.         // modulo is slow
  22.         if (ltid % (2*stride) == 0) {
  23.             smem[tid] += smem[ltid + stride];
  24.         }
  25.         __syncthreads();
  26.     }
  27.    
  28.     // write result for this block to global output memory
  29.     if (tid == 0) {
  30.         output[blockIdx.x] = smem[0];
  31.     }
  32. }
  33.  
  34. // consider that size = 128
  35. dim3 block_dim(64, 1, 1);
  36. dim3 grid_dim(2, 1, 1);
  37. // amount of shared memory used
  38. int smem_size = block.x * sizeof(int)
  39. // gpu_output
  40. cudaMalloc((void **)&gpu_output, grid.x*sizeof(int));
  41. // call kernel
  42. reduction_v1<<<grid_dim, block_dim, smem_size>>>(gpu_input, gpu_output, size);
  43.  
  44.  

Reste un dernier problème : comment faire la somme des éléments du tableau output ?

Plusieurs solutions peuvent être proposées :

6.4. Travail

Implantez les algorithmes de réduction en CUDA et tester la performance de ces algorithmes pour différentes tailles de tableaux :

__threadfence() is used to halt the current thread until all previous writes to shared and global memory are visible by other threads.

6.5. Bibliographie et sitographie