Part 4 : La réation de Gray Scott en GPU avec Cuda

À tout seigneur tout honeur, commençons par l'implémentation Cuda histoire d'avoir un premier point de comparaison avec la version CPU que nous avions abordé dans le dernier cours.

Je laisse le soin à mes collègues bien plus compétents que moi pour les implémentations en Sycl, Kokkos, Alpaka, OpenCL ou encore OpenACC.

Mais ce coup là, au lieu de commencer par l'implémentation (ce qu'il ne faut jamais faire, mais comme l'exemple de la section 3.3 est très simple ça ne nous a pas posé de problème), nous commencerons par regarder de plus près ce qui se trame dans un GPU pendant un calcul.

Chaque calcul est effectué sur la grille qui est elle-même composée de blocs, eux-mêmes composés de threads. La grille et les blocs qui la composent peuvent être à 1, 2 ou 3 dimensions. Rappelez-vous :

Device 0: "Quadro M2200"
	...
	Maximum number of threads per block:           1024
	Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
	Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
	Maximum memory pitch:                          2147483647 bytes

Pour simuler la réaction de Gray Scott, nous aurons besoin d'images (à deux dimensions donc). Il est donc mieux de répartir les calculs aussi sur deux dimensions. Les accès mémoire seront plus rapides car plus proches physiquement de la vraie position des données.

Le problème c'est que nous voulons calculer des images FULL-HD soit 1920x1080, ce qui est bien entendu plus grand que 1024, la taille maximal d'un bloc.

Dans l'exemple de la section 3.3, notre produit de Hadamard était à une dimension. Donc nous n'avions qu'à vérifier si le nombre d'éléments est plus grand que le nombre maximal de thread par bloc, si ce n'est pas le cas, on se contente de faire un bloc, sinon, on essaie de faire des blocs pleins :

	int dimGridX = 1;
	int dimBlockX = 1;
	if(nbElement > maxNbThreadPerBlockX){          //au lieu d'utiliser plusieurs warps d'un block, on utilise plusieurs block
		dimBlockX = maxNbThreadPerBlockX;
		dimGridX = nbElement/dimBlockX;
	}else{                                //sinon on utlise qu'un seul bloc
		dimBlockX = nbElement;
	}

Ce cas de figure ne fonctionne pas si le nombre d'éléments n'a pas le bon goût d'être un multiple du nombre de thread maximal par bloc.

Si on résume la situation, il va falloir préparer notre produit de Hadamard, car nous ne voulons pas oublier d'éléments dans notre simulation, et en plus, le faire à deux dimensions et non une ^{noteJ'espère que le lecteur perspicace n'a pas râler jusque là sinon il va se dire que ce cours est bancal.}.

Dans l'hypothèse où nous voulons des blocs carrés, ceux-ci ne peuvent pas dépasser $32 \times 32 = 1024$ sinon Cuda ne sera pas content.

Heureusement pour nous 1920 à le bon goût de se diviser par 32 : $\dfrac{1920}{32} = 60$ , pour 1080, ce n'est pas le cas, mais il a quand même le bon goût de se diviser par 30 : $\dfrac{1080}{30} = 36$ .

Si on résume la situation :

Nous aurons $36 \times 60 = 2160$ blocs, qui est bien plus petit que la taille maximale de la grille
Nous aurons $32 \times 30 = 960$ threads par blocs, qui est inférieur à (le maximum de threads par bloc) donc Cuda sera content. On remarque au passage que notre nombre de threads par bloc est proche du maximum de threads par bloc, donc nous utiliserons, a priori, plutôt bien le GPU. Et c'est un multiple de (par construction) donc l'utilisation des warps sera efficace.

Évidemment, cela a été fait exprès. Le fait qu'un GPU soit prévu pour faire de l'affichage sur un écran FULL-HD est une volonté téchnique, donc rien de surprenant là-dedans.

La figure 14 illustre le découpage en blocs d'une image Full-HD.

Figure 14 : Illustration du découpage en blocs des images Full-HD.

Note : Bien entendu, il est possible d'utiliser la même méthode que pour l'implémentation en C++, en changeant les calculs sur les bords de nos images, mais comme il s'agit de calcul sur GPU il est important de faire le plus possible de calculs identiques pour utiliser efficacement les blocs.

Nous allons donc utiliser un padding sur tous les bords de nos images pour que le calcul soit le même partout. Nous ferons donc nos calculs en 1922x1082 au lieu de 1920x1080, soit un "surcoup" de mémoire de $0.3\%$ qui paraît acceptable.

Il reste une dermière difficulté, la taille de la mémoire du GPU, dans mon cas environ 4 Go :

Device 0: "Quadro M2200"
	...
	Total amount of global memory:                 4044 MBytes (4240179200 bytes)

Or, nous voulons simuler $34\,000$ images et en conserver $1\,000$ (sachant que se sont des matrices de float). Donc $1\,920 \times 1\,080 \times 1\,000 = 2\,073\,600\,000$ float, $8\,294\,400\,000$ octets, soit environ 8.3

Go.

Si on s'arrête là, on pourrait se dire que cette simulation serait faite en 2 fois. Mais il faut prendre en compte la taille des temporaire que nous utiliserons.

Nous aurons besoin de 4 matrices temporaires ( U_1

) de taille 1922x1082 soit $1\,922 \times 1\,082 \times 4 = 8\,318\,416$ float soit $33\,273\,664$ octets.

Il reste la matrice $\nabla^2$ de

float, soit

octets :

Il faut reconnaitre que cette matrice n'est pas passionante, mais il est important d'être dans le même cas de calcul que dans le cours précédent où cette matrice existait.

Bref, en résumé :

Nous avons donc

Mo de marge sur notre simulation et nous pouvons donc faire ce calcul en 2 fois seulement. Ce calcul n'est valable que si aucun autre processus utilise le GPU (sur mon ordinateur, par exemple, je sais que l'OS utilise environ 315 Mo pour gérer l'affichage et deux trois autres bricoles du genre). L'API CUDA donne la mémoire totale, et non la mémoire disponible (jusqu'à ce qu'ils me fasse mentir avec une mise à jour, mais pour le moment ce n'est pas le cas).

Attention : un programme Cuda ne compile pas en C++17 avec nvcc (en tout cas avec la version 11.2). Et il ne vous donnera qu'une vieille erreur du genre :

/usr/include/c++/9/bits/stl_pair.h(442): error: argument list for class template "std::pair" is missing


/usr/include/c++/9/bits/stl_pair.h(442): error: expected a ")"


/usr/include/c++/9/bits/stl_pair.h(442): error: template parameter "_T1" may not be redeclared in this scope


/usr/include/c++/9/bits/stl_pair.h(442): error: expected a ";"

Donc faites attention. Il sait faire du C++11 mais pas du C++17.