Description¶

Dans le modèle de programmation Cuda, les threads sont réprésentés par une grille constituée de blocs de threads. La grille et les blocs sont de même dimension \(\mathtt{d}\) où \(\mathtt{d}\) est parmi 1, 2 ou 3. Le principe de tout ce qui va être expliqué ci-dessous est le même pour chacune des dimensions \(\mathtt{d}\). Puisque les deux exemples qui seront présentés dans ce document calculeront sur le GPU des coefficients d’une matrice, il sera plus commode de choisir \(\mathtt{d=2}\) mais \(\mathtt{d=1}\) ou \(\mathtt{d=3}\) pourrait s’adapter à nos exemples.

Dans le cas de la dimension 2, attention que l’ordre habituel \(\mathtt{(ligne, colonne)}\) dans lequel on se représente un tableau 2D se traduit par (ordonnée, abscisse) :

• une coordonnées en largeur est souvent notée en utilisant la lettre \(\mathtt{x}\) ; une largeur correspond à un écart de colonnes
• une coordonnées en hauteur est souvent notée en utilisant la lettre \(\mathtt{y}\) ; une hauteur correspond à un écart de lignes.

Entrons dans le détail :

• La grille est constituée de blocs à deux dimensions, \(\mathtt{nbx}\) blocs en largeur, \(\mathtt{nby}\) blocs en hauteur

• Chaque bloc de la grille est identifié par les deux indices

  • l’indice \(\mathtt{bx}\) de colonne dans la grille avec \(\mathtt{0\leq bx< nbx}\),
  • l’indice \(\mathtt{by}\) de ligne dans la grille avec \(\mathtt{0\leq by< nby}\).

• Chaque bloc est constitué de threads.

• Chaque bloc a deux dimensions \(\mathtt{ntx}\) et \(\mathtt{nty}\) et tous les blocs de la grille ont mêmes dimensions

• Chaque thread d’un bloc est référencé avec deux indices \(\mathtt{tx}\) et \(\mathtt{ty}\) où \(\mathtt{0\leq tx< ntx}\) et où \(\mathtt{0\leq ty< nty}\)

• Chaque bloc contient donc le même nombre de threads à savoir \(\mathtt{ntx \times nty}\).

• Le nombre de threads par bloc est matériellement limité, par exemple à 1024 ; on peut donc envisager des blocs de 32 x 32 threads.

• Si on voit la grille comme une grille de threads (et non de blocs) chaque thread a deux indices absolus :

  • l’indice de colonne : \(\mathtt{x = bx \times ntx + tx}\),
  • l’indice de ligne : \(\mathtt{y = by \times nty + ty}\).

Dans le dessin ci-dessus, on a une grille de \(\mathtt{nbx \times nby = 12 \times 4}\) blocs. Chaque bloc se décompose en \(\mathtt{ntx \times nty = 2 \times 4}\) threads. On a isolé, en bleu, dans le bloc (bx, by) = (8, 2) le thread (tx, ty) = (1, 2). La position absolue de ce thread est (x, y) = (17, 10).

Utlisation du modèle pour le calcul matriciel¶

Imaginons qu’on veuille effectuer le produit matriciel \(\mathtt{A\times B =: C}\) où

• \(\mathtt{A}\) est une matrice de taille lignes x colonnes \(\mathtt{(n, p)}\),
• \(\mathtt{B}\) est une matrice de taille lignes x colonnes \(\mathtt{(p, q)}\),
• \(\mathtt{C}\) est donc une matrice de taille lignes x colonnes \(\mathtt{(n, q)}\).

Pour fixer les idées, cf. figure ci-dessous où la matrice C est en bleu, on va supposer que

\(\mathtt{n=6}\), \(\mathtt{p=8}\), \(\mathtt{q=11}\)

en sorte que \(\mathtt{C}\) est de largeur 11 et de hauteur 6. Dans le modèle, on recouvre la matrice \(\mathtt{C}\) à calculer par une grille de blocs de threads. Pour fixer les idées, on va choisir chaque bloc de taille dimt= (2, 4) donc chaque bloc contient \(\mathtt{2\times 4=8}\) threads :

On posera \(\mathtt{ntx=2}\) et \(\mathtt{nty=4}\). Calculons le nombres de blocs en largeur et en hauteur pour recouvrir la matrice \(\mathtt{C}\) :

• chaque bloc à une largeur de \(\mathtt{ntx=2}\) et \(\mathtt{q=11}\) donc il faut en largeur \(\mathtt{nbx=6}\) blocs (et il y aura \(\mathtt{6\times 2 -11=1}\) thread inutile en largeur) ;
• chaque bloc à une hauteur de \(\mathtt{nty=4}\) et \(\mathtt{n=6}\) donc il faut en hauteur \(\mathtt{nby=2}\) blocs (et il y aura \(\mathtt{4\times 2 -6=2}\) threads inutiles en hauteur).

Comme on le voit sur la figure, la grille est formée de \(\mathtt{2 \times 6}\) blocs.

La grille ainsi constituée a une largeur de 12 threads et une hauteur de 8 threads et donc peut exécuter 96 threads, ce qui excède le nombre total des 66 coefficients de C à calculer.

Initialement, les matrices \(\mathtt{A}\) et \(\mathtt{B}\) ainsi qu’une matrice \(\mathtt{C}\) nulle et à la bonne taille sont définies sur l’hôte. Ces matrices sont envoyées sur le device. Pour calculer la matrice \(\mathtt{C}\), il y aura au moins autant de threads sur le GPU que de coefficients de \(\mathtt{C}\) à calculer. Dans notre exemple, il y aura 96 threads indexés \(\mathtt{(x, y)}\) avec \(\mathtt{0\leq x< 12}\) et \(\mathtt{0\leq y< 8}\).

Quand on lance un kernel, cf. les exemples de code plus bas, outre les données utiles au calcul, ici \(\mathtt{A}\), \(\mathtt{B}\) et \(\mathtt{C}\), on lui donne le couple \(\mathtt{dimb=(nbx, nby)}\) et le couple \(\mathtt{dimt=(ntx, nty)}\). Par ailleurs, quand ce kernel s’exécute, il sait dans quel bloc \(\mathtt{(bx, by)}\) de la grille il s’exécute et quel thread \(\mathtt{(tx, ty)}\) du bloc il exécute. Cela lui permet de savoir dans quel thread \(\mathtt{(x, y)}\) de la grille il s’exécute. Et notre kernel va être écrit en sorte que l’élément en colonne \(\mathtt{x}\) et ligne \(\mathtt{y}\) soit représenté par le thread indexé \(\mathtt{(x, y)}\).

En pratique, il est fréquent que les blocs de threads soient « carrés », ie que \(\mathtt{ntx=nty}\).

Références¶

Le modèle de programmation Cuda est décrit en détail dans un document officiel : Cuda programming model.

Plus généralement, concernant la programmation Cuda, on pourra consulter l’ouvrage progressif et détaillé CUDA C Programming paru en 2014.

Le kernel¶

Il est déclaré ligne 5. Il s’exécute sur le device. Il reçoit les matrices \(\mathtt{A}\), \(\mathtt{B}\) et \(\mathtt{C}\) qui se trouvent sur le device. La fonction est décorée ligne 4 par le décorateur cuda.jit. C’est ce décorateur qui indique à Numba que la fonction est un kernel. On passe à ce décorateur, sous forme de chaîne de caractères, ici

'(int64[:,:], int64[:,:], int64[:,:])'

un triplet de types de données transmises au kernel. Ici int64[:,:] désigne un tableau 2d d’entiers sur 64 bits. On trouvera quelques indications sur des signatures courantes ici : signature specifications.

Les kernels s’exécutent en parallèle. Chaque kernel calcule à la ligne 16 le coefficient de \(\mathtt{C}\) à la ligne \(\mathtt{y}\) et à la colonne \(\mathtt{x}\) de la grille. Ce calcul résulte tout simplement du produit ligne par colonne, effectué dans le kernel par la boucle for de la ligne 15.

Le code ci-dessus ne le montre pas, mais le kernel quand il s’exécute, connait le bloc et le thread d’exécution ; par ailleurs, suite à l’appel ligne 43, le kernel dispose des dimensions de chaque bloc (le dimt) et des dimensions de la grille en blocs (le dimb). Donc, à l’aide d’une formule assez simple, le kernel est en mesure de calculer la position absolue \(\mathtt{(x, y)}\) du thread dans la grille. Pour nous épargner le calcul, Numba met à disposition ligne 10 une fonction cuda.grid qui calcule cette position. Il faut donner à cette fonction la dimension des blocs, dans notre exemple, ce sont les 2 dimensions (le \(\mathtt{d=2}\) expliqué dans la description du modèle Cuda) mais on avait vu en introduction du modèle de programmation Cuda qu’il serait possible de choisir des dimensions 1 ou 3.

On a vu que la grille de threads déborde de la matrice \(\mathtt{C}\). Donc, le kernel doit d’assurer (cf. ligne 12) que le calcul effectué en \(\mathtt{(x, y)}\) correspond à une position légitime de la matrice \(\mathtt{C}\).

Un kernel ne renvoie rien, il agit par effet de bord.

Dimensions de la grille¶

Le nombre de threads qu’un bloc peut exécuter est plafonné. Par exemple, pour la carte graphique que j’ai utilisée, il est de 1024. Dans ce cas et si les blocs sont de dimension 2, il est courant de choisir un bloc de taille \(\mathtt{32\times 32}\). Dans notre petit exemple, on a choisi une taille \(\mathtt{(2, 4)}\), cf. lignes 56-57. Quand on connait les dimensions de chaque blocs, il est facile de calculer le nombre de blocs suivant chaque dimension pour recouvrir la matrice. C’est ce que fait la fonction dims(n, q, nt). L’idée du calcul par les formules lignes 30 et 33 a été donné lors de la description de l’exemple où on avait trouvé une grille de dimension de \(\mathtt{(2, 6)}\). Plus loin dans ce document, un aparté explique comment trouver cette formule et lui donner la forme habituellement rencontrée.

Calcul du produit matriciel¶

Les étapes sont les suivantes :

• lignes 19-20 et 21-22 : sur le host, on génère aléatoirement des matrices Numpy \(\mathtt{A}\) et \(\mathtt{B}\) à coefficients entre -10 et 10 et de type int64 ainsi qu’une matrice Numpy \(\mathtt{C}\) initialisée à 0
• lignes 39-41 : on déplace les 3 matrices sur le device ; Numba se charge pour nous de l’allocation mémoire et du transfert.
• ligne 43 : on lance le kernel en donnant tous les éléments utiles : les 3 matrices, bien sûr, et aussi, entre crochets, les dimensions de chaque bloc et de la grille.
• ligne 44 : la synchronisation va avoir pour effet que l’hôte va interrompre son fil d’exécution jusqu’à ce que tous les threads exécutant le kernel aient terminé.
• ligne 46 : copie la matrice depuis le device vers le host.
• lignes 63 et 65 : on peut alors accéder au résultat et par exemple comparer.

Variante d’écriture du kernel¶

Le kernel écrit dans le programme ci-dessus :

@cuda.jit('(int64[:,:], int64[:,:], int64[:,:])')
def cu_square_matrix_mul(A, B, C):

    n, p=A.shape
    _, q=B.shape

    x, y = cuda.grid(2)

    if x >= q or y >= n:
        return

    for i in range(p):
        C[y, x] += A[y, i] * B[i, x]

est parfois écrit de manière moins directe. En effet, Numba met à notre disposition la fonction cuda.grid mais il est possible de calculer explicitement \(\mathtt{x}\) et \(\mathtt{y}\). D’où le code suivant :

@cuda.jit('(int64[:,:], int64[:,:], int64[:,:])')
def cu_square_matrix_mul(A, B, C):

    n, p=A.shape
    _, q=B.shape

    tx = cuda.threadIdx.x
    ty = cuda.threadIdx.y
    bx = cuda.blockIdx.x
    by = cuda.blockIdx.y
    bw = cuda.blockDim.x
    bh = cuda.blockDim.y

    x = tx + bx * bw
    y = ty + by * bh

    if x >= q or y >= n:
        return

    for i in range(p):
        C[y, x] += A[y, i] * B[i, x]

En effet, chaque thread connait la position de son bloc dans la grille. On y a accède de la manière suivante :

bx = cuda.blockIdx.x
by = cuda.blockIdx.y

Il en est de même pour le thread dans son bloc :

tx = cuda.threadIdx.x
ty = cuda.threadIdx.y

Le thread connaît aussi la dimension nb de chaque bloc car il lui a été fournie entre crochets au moment du lancement du kernel :

cu_square_matrix_mul[nb, nt](dA, dB, dC)

Depuis le kernel, on y accède de la manière suivante :

bw = cuda.blockDim.x
bh = cuda.blockDim.y

autrement dit, on a nb = (bw, bh).

À partir de là, et vue la structure de la grille

on en déduit que

x = tx + bx * bw
y = ty + by * bh

Un autre algorithme ?¶

Pour information, le produit de deux matrices de taille 10000 calculé en un seul processus par un algorithme de Strassen optimisé met plus de 3 min 30 s. Voici le code sous SageMath :

# SageMath et lib flint2

n=10000
p=10000
q=10000

start=time()
A=random_matrix(ZZ,n,p, x=-10, y=10)
B=random_matrix(ZZ,p,q, x=-10, y=10)
delta = time() - start

print(f"Durée : {delta:.2f}s")



start=time()
A*B
delta = time() - start

print(f"Durée : {delta:.2f}s")

Durée : 18.86s
Durée : 219.95s

Alternative de démonstration¶

Ce qui suit n’a d’intérêt que du point de vue mathématique et peut-être ignoré. Il s’agit juste de montrer comment on pouvait trouver la formule directement.

Si \(\mathtt{n\leq d}\) alors \(\mathtt{N=1}\) . On suppose désormais que \(\mathtt{n>d}\) ce qui entraîne que \(\mathtt{N\geq 1}\). Alors comme \(\mathtt{N}\) est le plus petit entier vérifiant la propriété, \(\mathtt{N-1}\) ne la vérifie pas et comme \(\mathtt{N-1\geq 0}\), c’est que \(\mathtt{(N-1)d<n}\)

Montrons que \(\mathtt{N}\) est forcément le plus grand entier qui vérifie \(\mathtt{(N-1)d<n}\). En effet, puisque pour tout entier \(\mathtt{M>N-1}\) et donc, puisque \(\mathtt{M}\) et \(\mathtt{N}\) sont des entiers, forcément tels que \(\mathtt{M\geq N}\), on a \(\mathtt{Md\geq Nd\geq n}\).

Réciproquement, si \(\mathtt{N}\) est le plus grand entier qui vérifie \(\mathtt{(N-1)d<n}\) alors forcément, comme \(\mathtt{N+1}\) ne vérifie pas la condition, on a \(\mathtt{Nd\geq n}\) et \(\mathtt{N}\) est forcément le plus petit entier qui vérifie cette propriété puisque si \(\mathtt{M<N}\) alors (inégalité entre entiers) \(\mathtt{M\leq N-1}\) donc \(\mathtt{Md\leq (N-1)d<n}\).

Donc, finalement, on cherche le plus grand entier \(\mathtt{N}\) tel que \(\mathtt{(N-1)d<n}\) ou encore, puisque c’est une inégalité entre entiers, tel que \(\mathtt{(N-1)d\leq n-1}\) donc \(\mathtt{Nd\leq n+d-1}\) et on sait que le plus grand entier \(\mathtt{N}\) qui vérifie cela est le quotient de la division entière \(\mathtt{(n+d-1)//d}\) ce qui redonne le résultat vu plus haut.