cuBLAS et cuDNN

Pourquoi ne pas réécrire GEMM et convolution soi-même

Un ingénieur qui a suivi les chapitres précédents sait écrire un kernel de multiplication matricielle correct, avec tuiles en mémoire partagée et accès coalescés. C'est une étape pédagogique nécessaire, mais un kernel « fait main » plafonne en général entre 40 et 60 % du débit crête FP32 d'un GPUGPUIAProcesseur graphique parallélisant massivement les calculs matriciels ; indispensable à l'entraînement et à l'inférence des modèles de deep learning.Voir dans le glossaire récent, et beaucoup moins dès qu'on vise les Tensor Cores en FP16 ou TF32. cuBLAS et cuDNN exploitent des formes de tuiles et des séquences d'instructions MMA retunées à chaque génération de silicium ; réinventer cette couche a un coût d'opportunité élevé pour un gain presque toujours négatif.

Ce chapitre ne traite donc pas de « comment coder du BLAS », mais de ce qu'il faut savoir pour utiliser correctement ces deux bibliothèques : conventions de layout mémoire, éventail des algorithmes disponibles pour une même opération, méthode de sélection, et les rares situations où réécrire reste justifié.

cuBLAS hérite de la convention Fortran/BLAS historique : column-major. La quasi-totalité du code C++ et Python moderne (tableaux C classiques, NumPy par défaut, la plupart des frameworks) est row-major. Ignorer cette différence ne fait pas planter le programme : il produit un résultat silencieusement transposé, qui passe parfois inaperçu sur des cas de test mal choisis.

GEMM : le piège du layout column-major

Considérons un cas simple : multiplier deux matrices A (M×K) et B (K×N) stockées en row-major, comme le fait naturellement tout code C++ qui alloue float* A = new float[M*K] et remplit ligne par ligne. Appeler cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, M, N, K, &alpha, A, M, B, K, &beta, C, M) sans précaution ne produit pas A × B : cuBLAS interprète les mêmes octets comme des matrices column-major, donc comme si A était en réalité une matrice K×M et B une matrice N×K. Le résultat calculé n'est pas une erreur numérique, c'est un calcul parfaitement correct — mais sur des matrices différentes de celles que vous pensiez lui donner.

La solution standard ne nécessite aucune copie mémoire. Elle repose sur une identité algébrique simple :

(A × B)ᵀ = Bᵀ × Aᵀ

Une matrice row-major de dimensions M×N est, bit à bit, identique à sa transposée stockée en column-major de dimensions N×M — c'est le même buffer, seule l'interprétation change. On peut donc demander à cuBLAS de calculer Bᵀ × Aᵀ en lui passant directement les buffers B et A sans transposition explicite : le résultat, lu en column-major, est exactement C en row-major.

// A : M x K row-major, B : K x N row-major, C : M x N row-major
// on substitue A <-> B et M <-> N : cuBLAS calcule Ct = Bt x At,
// ce qui, relu en row-major, donne directement C = A x B
cublasSgemm(handle,
            CUBLAS_OP_N, CUBLAS_OP_N,
            N, M, K,
            &alpha,
            B, N,   // leading dimension de B : nombre de colonnes réelles
            A, K,   // leading dimension de A
            &beta,
            C, N);

Avec M=128, K=64, N=256, tout en row-major : l'appel correct inverse les rôles et passe N=256 puis M=128 en premiers paramètres de dimension, B avant A, et une leading dimension pour C égale à N=256 — la largeur réelle de C, pas M. Une équipe qui débogue un GEMM « qui donne des résultats bizarres uniquement sur des matrices rectangulaires » a, dans l'immense majorité des cas, oublié cette inversion ou mal renseigné une leading dimension.

La leading dimension (lda, ldb, ldc) est la deuxième source d'erreur classique, indépendante de la première. Elle ne représente pas la largeur logique de la matrice manipulée, mais le pas mémoire réel entre deux colonnes consécutives dans le buffer sous-jacent. Ce piège devient concret sur une sous-matrice extraite d'un tenseur plus large — un batch d'exemples, une matrice avec du padding d'alignement. Passer la largeur logique de la sous-matrice au lieu du pas réel du tableau parent ne provoque ni erreur de compilation ni plantage : cuBLAS lit et écrit à des offsets légèrement décalés, mélangeant silencieusement des donnéesdonnéesIAEnsemble d'informations structurées ou non utilisées pour entraîner, évaluer ou alimenter un modèle. La qualité, la quantité et la représentativité des données sont les facteurs décisifs pour les performances en apprentissage automatique.Voir dans le glossaire de lignes voisines.

Le paramètre beta de cublasSgemm contrôle l'accumulation : C = alpha * (A×B) + beta * C. Si beta n'est pas explicitement mis à 0 et que le buffer C n'a pas été préalablement initialisé (ou provient d'une itération précédente), le résultat inclut silencieusement les valeurs résiduelles de C. C'est une source fréquente de résultats numériques instables d'une exécution à l'autre, en particulier quand un buffer est réutilisé sans réinitialisation entre deux appels dans une boucle.

cuDNN : quatre façons de calculer une convolution

Une convolution 2D peut être exécutée par au moins quatre algorithmes distincts, avec des compromis mémoire/débit très différents. cuDNN les implémente tous et laisse le choix — explicite ou heuristique — à l'appelant.

im2col + GEMM. La méthode historique déplie chaque fenêtre de convolution en une colonne d'une matrice géante, ramenant le problème à une multiplication de matrices classique déléguée à cuBLAS. Simple conceptuellement, elle profite du débit crête GEMM, mais son coût mémoire explose : le buffer im2col peut peser plusieurs fois la taille du tenseur d'entrée, en particulier avec de grands noyaux ou un nombre élevé de canaux d'entrée.

Winograd. Cette famille d'algorithmes réduit le nombre de multiplications nécessaires au prix de davantage d'additions, via une transformation algébrique valable seulement pour des tailles de noyau spécifiques — typiquement 3×3, parfois 5×5. Le gain de débit est réel, jusqu'à environ 2× sur des noyaux 3×3 par rapport à un calcul direct. La contrepartie est une perte de précisionprécisionIAProportion des alertes émises par un modèle qui sont justifiées. Elle s'oppose au rappel : améliorer l'une dégrade l'autre.Voir dans le glossaire numérique légèrement supérieure, due aux transformations algébriques intermédiaires, un point à surveiller sur des réseaux sensibles à la reproductibilité bit-à-bit.

FFT et FFT-tiled. La convolution spatiale devient un produit terme-à-terme dans le domaine fréquentiel après transformation de Fourier. Le surcoût de la transformation elle-même rend cette approche rarement compétitive sur les noyaux 3×3 omniprésents en vision par ordinateur ; elle reprend l'avantage sur des noyaux plus larges, à partir de 7×7 environ.

Implicit GEMM. C'est l'algorithme par défaut sur les architectures récentes équipées de Tensor Cores (Volta et suivantes). Il reproduit l'effet d'im2col sans jamais matérialiser le buffer intermédiaire : chaque tuile est dépliée à la volée en mémoire partagée juste avant l'appel aux unités MMA. Il combine le débit d'un GEMM tuilé optimisé et une empreinte mémoire proche de celle d'un calcul direct — c'est pour cette raison qu'il domine la plupart des déploiements de production en inférenceinférenceIAUtilisation d'un modèle déjà entraîné sur une donnée nouvelle. Peu coûteuse à l'unité mais répétée à chaque requête, elle constitue le coût récurrent d'exploitation.Voir dans le glossaire comme en entraînement.

Comment cuDNN sélectionne un algorithme

cuDNN expose deux familles de fonctions pour choisir un algorithme. Les heuristiques (cudnnGetConvolutionForwardAlgorithm_v7) renvoient un classement basé sur des règles internes construites à partir de mesures antérieures, sans exécuter réellement le calcul sur le device : la réponse est quasi instantanée. La recherche exhaustive (cudnnFindConvolutionForwardAlgorithm) benchmarke chaque algorithme compatible directement sur le GPU cible et mesure le temps réel, ce qui donne un résultat plus fiable au prix d'un coût de démarrage — de quelques millisecondes à plusieurs secondes selon la taille du problème et le nombre de candidats testés.

Dans un pipeline avec des formes de tenseur fixes (batch size, résolution constants), le coût de la recherche exhaustive est amorti une seule fois puis mis en cache : c'est l'approche que PyTorch adopte via torch.backends.cudnn.benchmark = True. Si les formes varient à chaque appel — séquences de longueur variable, batchs de tailles hétérogènes — ce même réglage devient contre-productif : cuDNN relance un benchmark complet à chaque nouvelle forme rencontrée, ce qui peut ralentir l'ensemble du pipeline au lieu de l'accélérer.

Chaque algorithme a également un besoin en workspace — mémoire scratch temporaire allouée pour la durée de l'appel — qui varie fortement selon la méthode : implicit GEMM peut ne demander aucune mémoire supplémentaire, tandis qu'im2col peut réclamer plusieurs centaines de mégaoctets par couche pour un grand tenseur. cuDNN permet d'interroger ce besoin (cudnnGetConvolutionForwardWorkspaceSize) puis de le contraindre à un budget disponible ; si ce budget est trop serré, seuls les algorithmes économes en mémoire restent éligibles, parfois au prix d'un débit inférieur à ce que permettrait le matériel.

Quand réécrire un kernel plutôt qu'appeler la bibliothèque

Trois situations légitiment un kernel maison, et elles sont plus rares qu'on ne l'imagine en pratique.

1. Fusion d'opérateurs. Enchaîner convolution, ajout du biaisbiaisIARégularité correctement apprise dans des données qui ne représentent pas la réalité visée, ou qui enregistrent des décisions passées avec leurs préjugés. Changer d'algorithme ne le corrige pas.Voir dans le glossaire, activationfonction d'activationIAOpération non linéaire appliquée en sortie d'un neurone. Sans elle, empiler des couches serait inutile : une succession d'opérations linéaires reste équivalente à une seule.Voir dans le glossaire et normalisation dans un seul kernel évite plusieurs aller-retours en mémoire globale entre opérations séparées. cuDNN propose des API de fusion pour les motifs les plus courants (cudnnConvolutionBiasActivationForward), mais leur couverture reste limitée à un nombre restreint de combinaisons ; au-delà, un kernel personnalisé ou un compilateur de kernels comme Triton reprend la main.
2. Formes non standard. Convolutions groupées exotiques, noyaux creux, motifs d'accès irréguliers propres à un domaine spécifique — imagerie médicale volumique, traitement de graphes, signaux 1D à très haute fréquence d'échantillonnage. cuDNN couvre très bien les cas courants de la vision par ordinateur classique, nettement moins bien les cas de niche.
3. Exigence de précision ou de déterminisme strict. Certains algorithmes, Winograd en particulier, ne sont pas rigoureusement bit-reproductibles d'une exécution à l'autre sur le même matériel. En présence d'une contrainte de reproductibilité stricte — validation réglementaire, débogage d'un écart numérique entre deux runs — il faut forcer un algorithme déterministe, au prix d'un débit potentiellement inférieur.

En dehors de ces trois cas, réécrire une GEMM ou une convolution standard reste un pari structurellement perdant : le temps d'ingénierie investi dépasse largement le gain espéré, et chaque nouvelle génération de GPU imposerait de retuner le kernel maison, alors que cuBLAS et cuDNN se mettent à jour automatiquement avec la version du driver et de la toolkit.

Intégration côté C++

L'appel direct à l'API C de cuBLAS ou cuDNN nécessite trois objets vivant côté host mais liés au device. Un handle (cublasHandle_t, cudnnHandle_t) est créé une fois par thread ou par contextefenêtre de contexteIAQuantité de texte qu'un modèle peut prendre en compte simultanément : question, documents fournis et historique. Au-delà, les éléments les plus anciens sortent du champ.Voir dans le glossaire d'exécution — il porte l'état interne de la bibliothèque. Un stream CUDA est associé au handle via cublasSetStream ou cudnnSetStream, ce qui permet de recouvrir l'exécution avec d'autres opérations asynchrones plutôt que de tout sérialiser sur le stream par défaut. Des descripteurs décrivent enfin la forme exacte des données manipulées (cudnnTensorDescriptor_t, cudnnFilterDescriptor_t, cudnnConvolutionDescriptor_t), y compris le format mémoire (NCHW ou NHWC) — un mauvais choix ici peut à lui seul écarter les algorithmes les plus rapides, NHWC étant généralement préféré pour tirer parti des Tensor Cores.

Un cublasHandle_t ou cudnnHandle_t n'est pas thread-safe par défaut. Le partager entre plusieurs threads sans synchronisation externe explicite provoque des corruptions d'état intermittentes, difficiles à reproduire car dépendantes du timing d'exécution. Le pattern robuste consiste à créer un handle par thread (ou par flux logique de travail), pas un handle global partagé.

Intégration côté Python

En pratique, peu d'équipes appellent cuBLAS ou cuDNN directement depuis Python : PyTorch, TensorFlow ou cuPy encapsulent déjà ces bibliothèques et exposent les leviers utiles — choix d'algorithme, précision, budget de workspace — via des indicateurs haut niveau plutôt que des appels C bruts.

import torch

torch.backends.cudnn.benchmark = True        # recherche exhaustive, résultat mis en cache par forme
torch.backends.cuda.matmul.allow_tf32 = True # GEMM en TF32 sur Ampere et architectures suivantes
torch.backends.cudnn.deterministic = False   # autorise Winograd et les algorithmes non déterministes

Le recours direct à ctypes ou cffi pour appeler libcublas.so ou libcudnn.so depuis Python reste justifié dans deux cas précis : un besoin de contrôle fin absent des frameworks haut niveau (dimensionnement exact du workspace, stream personnalisé pour du multi-GPU asynchrone), ou un environnement de production sans framework MLapprentissage automatiqueIABranche de l'IA où le programme dégage lui-même ses règles à partir d'exemples, au lieu de les recevoir d'un développeur. C'est ce déplacement des règles écrites vers les régularités apprises qui définit l'IA moderne.Voir dans le glossaire complet — un service d'inférence minimal en C++ exposé via un binding Python léger avec pybind11.

Checklist avant mise en production

• Layout des données confirmé (row-major ou column-major) et vérifié sur un cas de test non symétrique — une matrice carrée symétrique masque souvent l'erreur.
• Leading dimensions validées séparément de la forme logique des matrices, en particulier sur des sous-matrices extraites d'un tenseur plus large.
• beta explicitement fixé à zéro pour tout appel GEMM où le buffer de sortie n'est pas censé accumuler une valeur préexistante.
• Algorithme cuDNN choisi en fonction du régime réel : shapes fixes → recherche exhaustive mise en cache ; shapes variables → heuristique.
• Budget de workspace dimensionné et vérifié en conditions de charge réelle, pas seulement en isolation sur une seule couche.
• Déterminisme évalué explicitement si le pipeline a une contrainte de reproductibilité (audit, débogage d'écart numérique).
• Format mémoire des tenseurs (NCHW vs NHWC) cohérent avec l'algorithme visé, en particulier pour tirer parti des Tensor Cores.

En résumé

cuBLAS et cuDNN ne sont pas de simples « boîtes noires rapides » : leur bon usage exige de comprendre la convention de layout attendue, l'éventail des algorithmes disponibles pour une opération donnée, et le compromis mémoire/débit/déterminisme que chaque choix implique. La compétence recherchée n'est pas de remplacer ces bibliothèques mais de savoir les piloter avec précision — et de reconnaître les rares cas où un kernel maison, ou une fusion via un compilateur comme Triton, reste réellement justifié.