Comment BLAS obtient-il des performances aussi extrêmes?

Par curiosité, j'ai décidé de comparer ma propre fonction de multiplication matricielle par rapport à l'implémentation BLAS ... J'ai été pour le moins surpris du résultat:

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

Took: 15.76542 seconds.

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

Took: 1.32432 seconds.

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

Ma mise en œuvre:

template<class ValT>
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("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

J'ai deux questions:

1. Étant donné qu'une multiplication matrice-matrice dit: nxm * mxn nécessite n * n * m multiplications, donc dans le cas au-dessus de 1000 ^ 3 ou 1e9 opérations. Comment est-il possible sur mon processeur 2,6 GHz que BLAS effectue des opérations 10 * 1e9 en 1,32 seconde? Même si les multiplcations étaient une seule opération et qu'il n'y avait rien d'autre à faire, cela devrait prendre ~ 4 secondes.
2. Pourquoi ma mise en œuvre est-elle tellement plus lente?

Un bon point de départ est le grand livre The Science of Programming Matrix Computations de Robert A. van de Geijn et Enrique S. Quintana-Ortí. Ils fournissent 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 n'opèrent que sur des vecteurs. Ces fonctions bénéficient de la vectorisation (par exemple de l'utilisation de SSE).

• Les fonctions de niveau 2 sont des opérations matrice-vecteur, par exemple un produit matrice-vecteur. Ces fonctions pourraient être implémentées en termes de fonctions de niveau 1. Cependant, vous pouvez améliorer les performances de ces fonctions si vous pouvez fournir une implémentation dédiée qui utilise 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 pouvez les implémenter en termes de fonctions Level2. Mais les fonctions Level3 effectuent des opérations O (N ^ 3) sur les 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 implémentation dédiée optimisée pour le cache / compatible avec le cache . Ceci est bien décrit dans le livre. Le principal avantage des fonctions Level3 provient de l'optimisation du cache. Cette augmentation dépasse largement la deuxième augmentation du parallélisme et d'autres optimisations matérielles.

À propos, la plupart (voire la totalité) des implémentations BLAS hautes performances ne sont PAS implémentées dans Fortran. ATLAS est implémenté en C. GotoBLAS / OpenBLAS est implémenté en C et ses pièces critiques de performance dans Assembler. Seule l'implémentation de référence de BLAS est implémentée dans Fortran. Cependant, toutes ces implémentations BLAS fournissent une interface Fortran telle qu'elle peut être liée à 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'a pas du tout d'importance).

IMHO no BLAS implémentation utilise des algorithmes comme l'algorithme Coppersmith – Winograd ou l'algorithme Strassen. Je ne suis pas exactement sûr de la raison, mais je suppose:

• Peut-être qu'il n'est pas possible de fournir une implémentation optimisée pour le cache de ces algorithmes (c'est-à-dire que vous perdriez plus que vous ne gagneriez)
• Ces algorithmes ne sont pas stables numériquement. Comme BLAS est le noyau de calcul de LAPACK, c'est un non-aller.

Modifier / mettre à jour:

Le nouveau document novateur sur ce sujet sont les papiers BLIS . Ils sont exceptionnellement bien écrits. Pour ma conférence "Bases du logiciel pour le calcul haute performance", j'ai implémenté le produit matrice-matrice en suivant leur article. En fait, j'ai implémenté plusieurs variantes du produit matrice-matrice. Les variantes les plus simples sont entièrement écrites en C brut et comportent moins de 450 lignes de code. Toutes les autres variantes optimisent simplement les boucles

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

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

ulmBLAS: Tutoriel sur GEMM (produit Matrix-Matrix)

Avec les papiers BLIS, il devient assez facile de comprendre comment des bibliothèques comme Intel MKL peuvent obtenir de telles performances. Et pourquoi peu importe que vous utilisiez le stockage principal en ligne ou en colonne!

Les benchmarks finaux sont ici (nous avons appelé notre projet ulmBLAS):

Benchmarks pour ulmBLAS, BLIS, MKL, openBLAS et Eigen

Une autre modification / mise à jour:

J'ai également écrit un tutoriel sur la façon dont BLAS est utilisé pour 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 d'étendre le tutoriel pour décrire et démontrer comment réaliser une implémentation parallèle hautement évolutive de la factorisation LU comme dans PLASMA .

Ok, c'est parti: Codage d'une factorisation LU parallèle optimisée pour le cache

PS: J'ai également fait des expériences pour améliorer les performances d'uBLAS. Il est en fait assez simple de booster (ouais, jouer sur les mots :)) les performances d'uBLAS:

Expériences sur uBLAS .

Voici un projet similaire avec BLAZE :

Expériences sur BLAZE .

Donc tout d'abord BLAS est juste une interface d'environ 50 fonctions. Il existe de nombreuses implémentations concurrentes de l'interface.

Tout d'abord, je mentionnerai des choses qui sont en grande partie sans rapport:

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

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

Ces opérations de petite dimension de taille fixe (appelées noyaux) sont codées en dur dans un code d'assemblage spécifique au processeur à l'aide de plusieurs fonctionnalités de processeur de leur cible:

• Instructions de style SIMD
• Parallélisme des niveaux d'instruction
• Conscience du cache

En outre, ces noyaux peuvent être exécutés en parallèle les uns par rapport aux autres en utilisant plusieurs threads (cœurs de processeur), dans le modèle de conception de réduction de carte typique.

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

(Astuce: c'est pourquoi si vous utilisez ATLAS, vous feriez mieux de créer et de régler la bibliothèque à la main pour votre machine particulière plutôt que d'utiliser une bibliothèque pré-construite.)

Premièrement, il existe des algorithmes plus efficaces pour la multiplication matricielle que celui que vous utilisez.

Deuxièmement, votre CPU peut faire 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 flottants ou 2 doubles. (bien sûr, ce chiffre n'est pas précis non plus, 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 utiliser des alias. Il existe des mots-clés ou des indicateurs spécifiques au compilateur que vous pouvez spécifier pour indiquer au compilateur qu'ils n'ont pas d'alias. Sinon, vous devez utiliser d'autres types que les pointeurs bruts, qui résolvent le problème.
• Vous écrasez le cache en effectuant un parcours naïf de chaque ligne / colonne des matrices d'entrée. Vous pouvez utiliser le blocage pour effectuer autant de travail que possible sur un bloc plus petit de la matrice, qui tient dans le cache du processeur, avant de passer au bloc suivant.
• Pour les tâches purement numériques, Fortran est quasiment imbattable, et C ++ demande beaucoup de persuasion pour atteindre une vitesse similaire. Cela peut être fait, et il y a quelques bibliothèques qui le démontrent (en utilisant généralement des modèles d'expression), mais ce n'est pas trivial, et cela n'arrive pas simplement .

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

La plupart des arguments à la deuxième question - assembleur, découpage en blocs, etc. (mais pas moins de N ^ 3 algorithmes, ils sont vraiment surdéveloppés) - jouent un rôle. Mais la faible vitesse de votre algorithme est essentiellement due à la taille de la matrice et à la disposition malheureuse des trois boucles imbriquées. Vos matrices sont si volumineuses qu'elles ne rentrent pas en même temps dans la mémoire cache. Vous pouvez réorganiser les boucles de manière à ce que le plus possible soit effectué sur une ligne du cache, ce qui réduit considérablement les rafraîchissements du cache (la division BTW en petits blocs a un effet analogique, mieux si les boucles sur les blocs sont disposées de la même manière). Une implémentation de modèle pour les matrices carrées suit. Sur mon ordinateur, sa consommation de temps était d'environ 1:10 par rapport à l'implémentation standard (comme la vôtre). En d'autres termes: ne programmez jamais une multiplication matricielle le long du "

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

Encore une remarque: cette implémentation est encore meilleure sur mon ordinateur que de tout remplacer par la routine BLAS cblas_dgemm (essayez-la sur votre ordinateur!). Mais beaucoup plus rapide (1: 4) appelle directement dgemm_ de la bibliothèque Fortran. Je pense que cette routine n'est en fait pas du Fortran mais du code assembleur (je ne sais pas ce qu'il y a dans la bibliothèque, je n'ai pas les sources). Je ne comprends pas du tout pourquoi cblas_dgemm n'est pas aussi rapide puisque, à ma connaissance, il s'agit simplement d'un wrapper pour dgemm_.

C'est une accélération réaliste. Pour un exemple de ce qui peut être fait avec l'assembleur SIMD sur le code C ++, voir quelques exemples de fonctions de matrice iPhone - elles étaient plus de 8x plus rapides que la version C, et ne sont même pas un assemblage "optimisé" - il n'y a pas encore de canalisation et là est des opérations de pile inutiles.

De plus, votre code n'est pas " restreint correct " - comment le compilateur sait-il que lorsqu'il modifie C, il ne modifie pas A et B?

En ce qui concerne le code d'origine dans MM multiply, la référence mémoire pour la plupart des opérations est la principale cause de mauvaises performances. La mémoire fonctionne 100 à 1000 fois plus lentement que le cache.

La plupart des accélérations proviennent de l'utilisation de techniques d'optimisation de boucle pour cette fonction triple boucle dans la multiplication MM. Deux techniques principales d'optimisation de boucle sont utilisées; déroulement et blocage. En ce qui concerne le déroulement, nous déroulons les deux boucles les plus externes et les bloquons pour la réutilisation des données dans le cache. Le déroulement de la boucle externe permet d'optimiser temporairement l'accès aux données en réduisant le nombre de références mémoire aux mêmes données à des moments différents pendant toute l'opération. Le blocage de l'index de boucle à un numéro spécifique aide à conserver les données dans le cache. Vous pouvez choisir d'optimiser pour le cache L2 ou le cache L3.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Pour de nombreuses raisons.

Premièrement, les compilateurs Fortran sont hautement optimisés et le langage leur permet de l'être. C et C ++ sont très lâches en termes de gestion des tableaux (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 obligé de créer du code générique. Dans Fortran, vos cas sont plus rationalisés, et le compilateur a un meilleur contrôle de ce qui se passe, ce qui lui permet d'optimiser davantage (par exemple en utilisant des registres).

Une autre chose est que Fortran stocke les éléments par colonne, tandis que C stocke les données par ligne. Je n'ai pas vérifié votre code, mais faites attention à la façon dont vous exécutez le produit. En C, vous devez scanner les lignes: de cette façon, vous scannez votre tableau le long de la mémoire contiguë, réduisant ainsi les échecs de cache. L'absence 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 dans l'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, la plupart répétées 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 des opérations coûteuses.

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

template<class ValT>
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("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

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

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