Comment BLAS parvient-il à une telle performance extrême ?

Question

Par curiosité, j'ai décidé de comparer les performances de ma fonction de multiplication matricielle personnelle avec l'implémentation BLAS... J'ai été pour le moins surpris du résultat :

Implémentation personnalisée, 10 essais de multiplication matricielle de 1000x1000 :

A pris : 15.76542 secondes.

Implémentation BLAS, 10 essais de multiplication matricielle de 1000x1000 :

A pris : 1.32432 secondes.

Ceci utilise des nombres à virgule flottante simple précision.

Ma Implémentation :

template
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Erreur tailles incorrectes");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2

`

J'ai deux questions :

1. Étant donné qu'une multiplication matricielle, par exemple : nxm mxn nécessite nnm multiplications, donc dans le cas ci-dessus 1000^3 ou 1e9 opérations. Comment est-il possible sur mon processeur 2.6Ghz pour BLAS de faire 101e9 opérations en 1.32 secondes ? Même si les multiplications étaient une seule opération et qu'il n'y avait rien d'autre à faire, cela devrait prendre environ ~4 secondes.
2. Pourquoi mon implémentation est-elle beaucoup plus lente ?

`

Answer 1

Un bon point de départ est le excellent livre The Science of Programming Matrix Computations par Robert A. van de Geijn et Enrique S. Quintana-Ortí. Ils proposent une version téléchargeable gratuitement.

BLAS est divisé en trois niveaux :

• Le niveau 1 définit un ensemble de fonctions d'algèbre linéaire qui opèrent uniquement sur des vecteurs. Ces fonctions bénéficient de la vectorisation (par exemple à partir de l'utilisation de SIMD tels que SSE).

• Les fonctions du niveau 2 sont des opérations matrice-vecteur, par exemple un produit matrice-vecteur. Ces fonctions pourraient être mises en œuvre en termes de fonctions de niveau 1. Cependant, vous pouvez améliorer les performances de ces fonctions si vous pouvez fournir une mise en œuvre dédiée qui exploite une architecture multiprocesseur avec mémoire partagée.

• Les fonctions de niveau 3 sont des opérations comme le produit matrice-matrice. Encore une fois, vous pourriez les implémenter en termes de fonctions de niveau 2. Mais les fonctions de niveau 3 effectuent des opérations O(N^3) sur des données O(N^2). Donc, si votre plate-forme dispose d'une hiérarchie de cache, vous pouvez améliorer les performances si vous fournissez une mise en œuvre dédiée qui est optimisée pour le cache/amie du cache. Cela est bien décrit dans le livre. L'augmentation principale des fonctions de niveau 3 provient de l'optimisation du cache. Cette amélioration dépasse de loin la deuxième augmentation obtenue grâce à la parallélisme et aux autres optimisations matérielles.

Soit dit en passant, la plupart (voire toutes) des implémentations BLAS haute performance ne sont PAS faites en Fortran. ATLAS est implémenté en C. GotoBLAS/OpenBLAS est implémenté en C et ses parties critiques en termes de performances en assembleur. Seule la mise en œuvre de référence de BLAS est faite en Fortran. Cependant, toutes ces implémentations BLAS fournissent une interface Fortran permettant de lier avec LAPACK (LAPACK tire toutes ses performances de BLAS).

Les compilateurs optimisés jouent un rôle mineur à cet égard (et pour GotoBLAS/OpenBLAS, le compilateur n'est pas du tout important).

À mon avis, aucune implémentation BLAS n'utilise des algorithmes comme l'algorithme de Coppersmith–Winograd ou l'algorithme de Strassen. Les raisons probables sont :

• Peut-être qu'il n'est pas possible de fournir une mise en œuvre optimisée pour le cache de ces algorithmes (c'est-à-dire que vous perdriez plus que vous ne gagneriez)
• Ces algorithmes ne sont pas numériquement stables. Comme BLAS est le noyau computationnel de LAPACK, cela ne peut pas être accepté.
• Bien que ces algorithmes aient une belle complexité temporelle sur papier, la notation Big O cache une grande constante, donc cela ne devient viable que pour des matrices extrêmement grandes.

Modifier/Mettre à jour :

Les nouveaux et révolutionnaires documents sur ce sujet sont les documents BLIS. Ils sont exceptionnellement bien écrits. Pour mon cours "Bases Logicielles pour le Calcul Haute Performance", j'ai implémenté le produit matrice-matrice suivant leur document. En fait, j'ai implémenté plusieurs variantes du produit matrice-matrice. La variante la plus simple est entièrement écrite en C pur et comporte moins de 450 lignes de code. Toutes les autres variantes ne font qu'optimiser les boucles

    for (l=0; l

`

Les performances globales du produit matrice-matrice dépendent seulement de ces boucles. Environ 99,9% du temps est passé ici. Dans les autres variantes, j'ai utilisé des fonctions intrinsèques et du code assembleur pour améliorer les performances. Vous pouvez voir le tutoriel sur toutes les variantes ici :

ulmBLAS : Tutoriel sur GEMM (Produit Matrice-Matrice)

Avec les documents BLIS, il devient assez facile de comprendre comment des bibliothèques comme Intel MKL peuvent atteindre de telles performances. Et pourquoi cela n'a pas d'importance si vous utilisez un stockage par ligne ou par colonne !

Les derniers benchmarks se trouvent ici (nous avons appelé notre projet ulmBLAS) :

Benchmarks pour ulmBLAS, BLIS, MKL, openBLAS et Eigen

Autre modification/Mise à jour :

J'ai également écrit quelques tutoriels sur la façon dont BLAS est utilisée pour résoudre des problèmes d'algèbre linéaire numérique comme la résolution d'un système d'équations linéaires :

Factorisation LU à Haute Performance

(Cette factorisation LU est par exemple utilisée par Matlab pour résoudre un système d'équations linéaires.)

J'espère trouver le temps pour étendre le tutoriel afin de décrire et de montrer comment réaliser une implémentation parallèle hautement évolutive de la factorisation LU comme dans PLASMA.

D'accord, voici : Coder une Factorisation LU Parallèle Optimisée pour le Cache

P.S. : J'ai également réalisé des expériences pour améliorer les performances de uBLAS. En fait, il est assez simple d'améliorer (oui, jeu de mots :) ) les performances de uBLAS :

Expériences sur uBLAS.

Voici un projet similaire avec BLAZE :

Expériences sur BLAZE.

`

Answer 2

Tout d'abord, BLAS est simplement une interface d'environ 50 fonctions. Il existe de nombreuses implémentations concurrentes de l'interface.

Tout d'abord, je mentionnerai des choses largement non liées :

• Fortran vs C, cela ne fait aucune différence
• Algorithmes matriciels avancés tels que Strassen, les implémentations ne les utilisent pas car ils n'aident pas en pratique

La plupart des implémentations divisent chaque opération en petites opérations de matrices ou de vecteurs de petite dimension de manière plus ou moins évidente. Par exemple, une grande multiplication de matrices de 1000x1000 peut être divisée en une séquence de multiplications de matrices de 50x50.

Ces opérations de petite dimension de taille fixe (appelées noyaux) sont codées en dur dans du code d'assemblage spécifique au CPU utilisant plusieurs caractéristiques du CPU cible :

• Instructions de style SIMD
• Parallélisme au niveau de l'instruction
• Sensibilisation au cache

De plus, ces noyaux peuvent être exécutés en parallèle les uns par rapport aux autres en utilisant plusieurs threads (cœurs CPU), dans le schéma de conception typique de map-reduce.

Jetez un œil à ATLAS qui est l'implémentation BLAS open source la plus couramment utilisée. Il existe de nombreux noyaux concurrents différents, et pendant le processus de construction de la bibliothèque ATLAS, il lance une compétition entre eux (certains sont même paramétrés, de sorte que le même noyau peut avoir des paramètres différents). Il essaye différentes configurations puis sélectionne la meilleure pour le système cible particulier.

(Conseil : c'est pourquoi si vous utilisez ATLAS, il vaut mieux construire et ajuster la bibliothèque manuellement pour votre machine particulière plutôt que d'utiliser une précompilée.)

Answer 3

Tout d'abord, il existe des algorithmes plus efficaces pour la multiplication de matrices que celui que vous utilisez.

Ensuite, votre CPU peut effectuer beaucoup plus d'une instruction à la fois.

Votre CPU exécute 3-4 instructions par cycle, et si les unités SIMD sont utilisées, chaque instruction traite 4 nombres flottants ou 2 doubles. (bien sûr, ce chiffre n'est pas non plus précis, car le CPU ne peut généralement traiter qu'une seule instruction SIMD par cycle)

Troisièmement, votre code est loin d'être optimal :

• Vous utilisez des pointeurs bruts, ce qui signifie que le compilateur doit supposer qu'ils peuvent être connectés. Il existe des mots-clés ou des indicateurs spécifiques au compilateur que vous pouvez spécifier pour dire au compilateur qu'ils ne sont pas connectés. Alternativement, vous devriez utiliser d'autres types que des pointeurs bruts, qui résolvent le problème.
• Vous épuisez le cache en effectuant une traversée naïve de chaque ligne/colonne des matrices d'entrée. Vous pouvez utiliser le blocking pour effectuer autant de travail que possible sur un bloc plus petit de la matrice, qui tient dans le cache du CPU, avant de passer au bloc suivant.
• Pour les tâches purement numériques, Fortran est pratiquement imbattable, et C++ demande beaucoup d'efforts pour atteindre une vitesse similaire. Cela peut être fait, et il y a quelques bibliothèques le démontrant (généralement en utilisant des modèles d'expressions), mais ce n'est pas trivial, et cela ne se passe pas simplement.

Answer 4

Je ne connais pas spécifiquement l’implémentation BLAS mais il existe des algorithmes plus efficaces pour la multiplication de matrices qui ont une meilleure complexité que O(n3). Un bien connu est l'algorithme de Strassen

Answer 5

Pour de nombreuses raisons.

Tout d'abord, les compilateurs Fortran sont hautement optimisés, et le langage leur permet d'être ainsi. C et C++ sont très libres en termes de manipulation de tableau (par exemple, le cas des pointeurs se référant à la même zone mémoire). Cela signifie que le compilateur ne peut pas savoir à l'avance quoi faire, et est contraint de créer du code générique. En Fortran, vos cas sont plus rationalisés, et le compilateur a un meilleur contrôle de ce qui se passe, lui permettant d'optimiser davantage (par exemple, en utilisant des registres).

Autre chose est que Fortran stocke les données colonne par colonne, tandis que C stocke les données ligne par ligne. Je n'ai pas vérifié votre code, mais faites attention à la façon dont vous effectuez le produit. En C, vous devez balayer ligne par ligne : de cette façon, vous balayez votre tableau le long d'une mémoire contiguë, réduisant les erreurs de cache. L'erreur de cache est la première source d'inefficacité.

Troisièmement, cela dépend de l'implémentation blas que vous utilisez. Certaines implémentations peuvent être écrites en assembleur et optimisées pour le processeur spécifique que vous utilisez. La version netlib est écrite en fortran 77.

De plus, vous effectuez de nombreuses opérations, dont la plupart sont répétitives et redondantes. Toutes ces multiplications pour obtenir l'indice sont préjudiciables à la performance. Je ne sais pas vraiment comment cela se fait dans BLAS, mais il existe de nombreuses astuces pour éviter les opérations coûteuses.

Par exemple, vous pourriez retravailler votre code de cette façon

template
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Erreur de tailles");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2

`

Essayez-le, je suis sûr que vous sauverez quelque chose.

Sur votre question n°1, la raison en est que la multiplication matricielle échelle en O(n^3) si vous utilisez un algorithme trivial. Il existe des algorithmes qui échellent beaucoup mieux.

`