Sous-programmes d'algèbre linéaire de base

Les sous-programmes d'algèbre linéaire de base ( BLAS ) sont une spécification qui prescrit un ensemble de routines de bas niveau pour effectuer des opérations d'algèbre linéaire courantes telles que l'addition de vecteurs , la multiplication scalaire , les produits scalaires , les combinaisons linéaires et la multiplication de matrices . Ce sont les routines de bas niveau standard de facto pour les bibliothèques d'algèbre linéaire ; les routines ont des liaisons pour C (« interface CBLAS ») et Fortran (« interface BLAS »). Bien que la spécification BLAS soit générale, les implémentations BLAS sont souvent optimisées pour la vitesse sur une machine particulière, donc leur utilisation peut apporter des avantages substantiels en termes de performances. Les implémentations BLAS tireront parti du matériel à virgule flottante spécial tel que les registres vectoriels ou les instructions SIMD .

Il s'agit d'une bibliothèque Fortran créée en 1979 et son interface a été normalisée par le Forum technique BLAS (BLAST), dont le dernier rapport BLAS est disponible sur le site Web de Netlib . Cette bibliothèque Fortran est connue comme l' implémentation de référence (parfois appelée de manière déroutante la bibliothèque BLAS) et n'est pas optimisée pour la vitesse mais est dans le domaine public .

La plupart des bibliothèques qui proposent des routines d'algèbre linéaire sont conformes à l'interface BLAS, permettant aux utilisateurs de la bibliothèque de développer des programmes indépendants de la bibliothèque BLAS utilisée.

De nombreuses bibliothèques BLAS ont été développées, ciblant différentes plateformes matérielles. Parmi les exemples, citons cuBLAS (NVIDIA GPU, GPGPU ), rocBLAS (AMD GPU) et OpenBLAS . Parmi les exemples de branches de bibliothèques BLAS basées sur le processeur, on trouve : OpenBLAS , BLIS (BLAS-like Library Instantiation Software) , Arm Performance Libraries, ATLAS et Intel Math Kernel Library (iMKL). AMD maintient un fork de BLIS optimisé pour la plateforme AMD . ATLAS est une bibliothèque portable qui s'optimise automatiquement pour une architecture arbitraire. iMKL est une bibliothèque gratuite et propriétaire optimisée pour x86 et x86-64 avec un accent sur les performances des processeurs Intel . OpenBLAS est une bibliothèque open source optimisée manuellement pour de nombreuses architectures populaires. Les tests de performance LINPACK s'appuient fortement sur la routine BLAS gemmpour ses mesures de performances.

De nombreuses applications logicielles numériques utilisent des bibliothèques compatibles BLAS pour effectuer des calculs d'algèbre linéaire, notamment LAPACK , LINPACK , Armadillo , GNU Octave , Mathematica , MATLAB , NumPy , R , Julia et Lisp-Stat.

Arrière-plan

Avec l'avènement de la programmation numérique, des bibliothèques de sous-routines sophistiquées sont devenues utiles. Ces bibliothèques contenaient des sous-routines pour des opérations mathématiques courantes de haut niveau telles que la recherche de racines, l'inversion de matrices et la résolution de systèmes d'équations. Le langage de choix était le FORTRAN . La bibliothèque de programmation numérique la plus importante était le Scientific Subroutine Package (SSP) d' IBM . Ces bibliothèques de sous-routines permettaient aux programmeurs de se concentrer sur leurs problèmes spécifiques et d'éviter de réimplémenter des algorithmes bien connus. Les routines de la bibliothèque étaient également meilleures que les implémentations moyennes ; les algorithmes matriciels, par exemple, pouvaient utiliser le pivotement complet pour obtenir une meilleure précision numérique. Les routines de la bibliothèque étaient également plus efficaces. Par exemple, une bibliothèque peut inclure un programme pour résoudre une matrice triangulaire supérieure. Les bibliothèques incluaient des versions simple précision et double précision de certains algorithmes.

Au départ, ces sous-routines utilisaient des boucles codées en dur pour leurs opérations de bas niveau. Par exemple, si une sous-routine devait effectuer une multiplication de matrice, elle aurait alors trois boucles imbriquées. Les programmes d'algèbre linéaire ont de nombreuses opérations de bas niveau courantes (les opérations dites « de noyau », sans rapport avec les systèmes d'exploitation ). Entre 1973 et 1977, plusieurs de ces opérations de noyau ont été identifiées. Ces opérations de noyau sont devenues des sous-routines définies que les bibliothèques mathématiques pouvaient appeler. Les appels de noyau présentaient des avantages par rapport aux boucles codées en dur : la routine de la bibliothèque était plus lisible, il y avait moins de risques de bugs et l'implémentation du noyau pouvait être optimisée pour la vitesse. Une spécification pour ces opérations de noyau utilisant des scalaires et des vecteurs , les Basic Linear Algebra Subroutines (BLAS) de niveau 1, a été publiée en 1979. BLAS a été utilisé pour implémenter la bibliothèque de sous-routines d'algèbre linéaire LINPACK .

L'abstraction BLAS permet une personnalisation pour des performances élevées. Par exemple, LINPACK est une bibliothèque à usage général qui peut être utilisée sur de nombreuses machines différentes sans modification. LINPACK pourrait utiliser une version générique de BLAS. Pour gagner en performances, différentes machines pourraient utiliser des versions personnalisées de BLAS. À mesure que les architectures informatiques sont devenues plus sophistiquées, des machines vectorielles sont apparues. BLAS pour une machine vectorielle pouvait utiliser les opérations vectorielles rapides de la machine. (Alors que les processeurs vectoriels ont fini par tomber en désuétude, les instructions vectorielles des processeurs modernes sont essentielles pour des performances optimales dans les routines BLAS.)

D'autres fonctionnalités de la machine sont devenues disponibles et pouvaient également être exploitées. Par conséquent, BLAS a été augmenté de 1984 à 1986 avec des opérations de noyau de niveau 2 qui concernaient les opérations vecteur-matrice. La hiérarchie de la mémoire a également été reconnue comme quelque chose à exploiter. De nombreux ordinateurs ont une mémoire cache beaucoup plus rapide que la mémoire principale ; garder les manipulations de matrice localisées permet une meilleure utilisation du cache. En 1987 et 1988, les BLAS de niveau 3 ont été identifiés pour effectuer des opérations matrice-matrice. Le BLAS de niveau 3 a encouragé les algorithmes partitionnés en blocs. La bibliothèque LAPACK utilise le BLAS de niveau 3.

Le BLAS original ne concernait que les vecteurs et les matrices densément stockés. D'autres extensions du BLAS, par exemple pour les matrices creuses, ont été abordées.

Fonctionnalité

Les fonctionnalités BLAS sont classées en trois ensembles de routines appelées « niveaux », qui correspondent à la fois à l'ordre chronologique de définition et de publication, ainsi qu'au degré du polynôme dans la complexité des algorithmes ; les opérations BLAS de niveau 1 prennent généralement un temps linéaire , O ( n ) , les opérations de niveau 2 un temps quadratique et les opérations de niveau 3 un temps cubique. Les implémentations BLAS modernes fournissent généralement les trois niveaux.

Niveau 1

Ce niveau comprend toutes les routines décrites dans la présentation originale de BLAS (1979), qui définissait uniquement des opérations vectorielles sur des tableaux à pas : produits scalaires , normes vectorielles , une addition vectorielle généralisée de la forme

{\displaystyle {\boldsymbol {y}}\leftarrow \alpha {\boldsymbol {x}}+{\boldsymbol {y <semantics> <mrow> <mstyle> <mrow> <mi>et</mi> </mrow> <mo>←</mo> <mi>α</mi> <mrow> <mi>x</mi> </mrow> <mo>+</mo> <mrow> <mi>et</mi> </mrow> </mstyle> </mrow> <annotation>{ displaystyle { boldsymbol {y}} leftarrow alpha { boldsymbol {x}} + { boldsymbol {y}}}</annotation> </semantics> </math></span><img src=

(appelé " axpy", "ax plus y") et plusieurs autres opérations.

Niveau 2

Ce niveau contient des opérations matrice-vecteur incluant, entre autres, une multiplication matrice - vecteur généralisée ( ) :gemv

{\displaystyle {\boldsymbol {y}}\leftarrow \alpha {\boldsymbol {A}}{\boldsymbol {x}}+\beta {\boldsymbol {y <semantics> <mrow> <mstyle> <mrow> <mi>et</mi> </mrow> <mo>←</mo> <mi>α</mi> <mrow> <mi>UN</mi> </mrow> <mrow> <mi>x</mi> </mrow> <mo>+</mo> <mi>β</mi> <mrow> <mi>et</mi> </mrow> </mstyle> </mrow> <annotation>{\displaystyle {\boldsymbol {y}}\leftarrow \alpha {\boldsymbol {A}}{\boldsymbol {x}}+\beta {\boldsymbol {y}}}</annotation> </semantics> </math></span><img src=

ainsi qu'un solveur pour x dans l'équation linéaire

{\displaystyle {\boldsymbol {T}}{\boldsymbol {x}}={\boldsymbol {y <semantics> <mrow> <mstyle> <mrow> <mi>T</mi> </mrow> <mrow> <mi>x</mi> </mrow> <mo>=</mo> <mrow> <mi>et</mi> </mrow> </mstyle> </mrow> <annotation>{\displaystyle {\boldsymbol {T}}{\boldsymbol {x}}={\boldsymbol {y}}}</annotation> </semantics> </math></span><img src=

avec T étant triangulaire. La conception du BLAS de niveau 2 a commencé en 1984, avec des résultats publiés en 1988. Les sous-routines de niveau 2 sont spécialement destinées à améliorer les performances des programmes utilisant BLAS sur des processeurs vectoriels , où les BLAS de niveau 1 sont sous-optimaux « car ils cachent la nature matricielle-vectorielle des opérations au compilateur. »

Niveau 3

Ce niveau, publié officiellement en 1990, contient des opérations matrice-matrice , y compris une « multiplication matricielle générale » ( gemm), de la forme

{\displaystyle {\boldsymbol {C}}\leftarrow \alpha {\boldsymbol {A}}{\boldsymbol {B}}+\beta {\boldsymbol {C}}, <semantics> <mrow> <mstyle> <mrow> <mi>C</mi> </mrow> <mo>←</mo> <mi>α</mi> <mrow> <mi>UN</mi> </mrow> <mrow> <mi>B</mi> </mrow> <mo>+</mo> <mi>β</mi> <mrow> <mi>C</mi> </mrow> <mo>,</mo> </mstyle> </mrow> <annotation>{\displaystyle {\boldsymbol {C}}\leftarrow \alpha {\boldsymbol {A}}{\boldsymbol {B}}+\beta {\boldsymbol {C}},}</annotation> </semantics> </math></span><img src=

où A et B peuvent éventuellement être transposés ou conjugués hermitiens à l'intérieur de la routine, et les trois matrices peuvent être parcourues. La multiplication matricielle ordinaire AB peut être effectuée en définissant α à un et C à une matrice entièrement composée de zéros de taille appropriée.

Le niveau 3 comprend également des routines de calcul

{\displaystyle {\boldsymbol {B}}\leftarrow \alpha {\boldsymbol {T}}^{-1}{\boldsymbol {B}}, <semantics> <mrow> <mstyle> <mrow> <mi>B</mi> </mrow> <mo>←</mo> <mi>α</mi> <msup> <mrow> <mi>T</mi> </mrow> <mrow> <mo>−</mo> <mn>1</mn> </mrow> </msup> <mrow> <mi>B</mi> </mrow> <mo>,</mo> </mstyle> </mrow> <annotation>{\displaystyle {\boldsymbol {B}}\leftarrow \alpha {\boldsymbol {T}}^{-1}{\boldsymbol {B}},}</annotation> </semantics> </math></span><img src=

où T est une matrice triangulaire , entre autres fonctionnalités.

En raison de l'omniprésence des multiplications de matrices dans de nombreuses applications scientifiques, y compris pour la mise en œuvre du reste du BLAS de niveau 3, et parce qu'il existe des algorithmes plus rapides au-delà de la répétition évidente de la multiplication matrice-vecteur, gemmc'est une cible de choix pour l'optimisation des implémenteurs de BLAS. Par exemple, en décomposant l'un ou les deux de A , B en matrices de blocs , gemmon peut implémenter de manière récursive . C'est l'une des motivations pour inclure le paramètre β , afin que les résultats des blocs précédents puissent être accumulés. Notez que cette décomposition nécessite le cas particulier β = 1 pour lequel de nombreuses implémentations optimisent, éliminant ainsi une multiplication pour chaque valeur de C. Cette décomposition permet une meilleure localité de référence à la fois dans l'espace et dans le temps des données utilisées dans le produit. Cela, à son tour, tire parti du cache sur le système. Pour les systèmes avec plus d'un niveau de cache, le blocage peut être appliqué une deuxième fois à l'ordre dans lequel les blocs sont utilisés dans le calcul. Ces deux niveaux d'optimisation sont utilisés dans des implémentations telles qu'ATLAS . Plus récemment, les implémentations de Kazushige Goto ont montré que le blocage uniquement pour le cache L2 , combiné à un amortissement prudent de la copie vers la mémoire contiguë pour réduire les échecs TLB , est supérieur à ATLAS . Une implémentation hautement optimisée basée sur ces idées fait partie de GotoBLAS , OpenBLAS et BLIS .

Une variante courante de gemmest le gemm3m, qui calcule un produit complexe en utilisant « trois multiplications de matrices réelles et cinq additions de matrices réelles au lieu des quatre multiplications de matrices réelles et deux additions de matrices réelles conventionnelles », un algorithme similaire à l'algorithme de Strassen décrit pour la première fois par Peter Ungar.

Implémentations

Accélérer: Le framework d' Apple pour macOS et iOS , qui inclut des versions optimisées de BLAS et LAPACK .
Bibliothèques de performances Arm: Bibliothèques de performances Arm, prenant en charge les processeurs Arm 64 bits basés sur AArch64 , disponibles sur Arm .
ATLAS: Logiciel d'algèbre linéaire réglé automatiquement , une implémentation open source des API BLAS pour C et Fortran 77. [
BLIS: Logiciel d'instanciation de bibliothèque de type BLAS pour une instanciation rapide. Optimisé pour la plupart des processeurs modernes. BLIS est une refactorisation complète de GotoBLAS qui réduit la quantité de code à écrire pour une plateforme donnée.
C++ et BLAS: La bibliothèque C++ AMP BLAS est une implémentation open source de BLAS pour l'extension de langage AMP de Microsoft pour Visual C++.
cuBLAS: BLAS optimisé pour les cartes GPU basées sur NVIDIA, nécessitant quelques appels de bibliothèque supplémentaires.
NVBLAS: BLAS optimisé pour les cartes GPU basées sur NVIDIA, fournissant uniquement des fonctions de niveau 3, mais en remplacement direct des autres bibliothèques BLAS.
clBLAS: Une implémentation OpenCL de BLAS par AMD. Fait partie des bibliothèques de calcul AMD.
clBLAST: Une implémentation OpenCL optimisée de la plupart des API BLAS.
BLAS propre: Une bibliothèque Fortran 77 et C BLAS implémentée sur la bibliothèque Eigen sous licence MPL , prenant en charge les architectures x86 , x86-64 , ARM (NEON) et PowerPC .
ESSL: Bibliothèque de sous-routines d'ingénierie et scientifiques d' IBM , prenant en charge l'architecture PowerPC sous AIX et Linux .
Aller à BLAS: Implémentation de BLAS sous licence BSD de Kazushige Goto , optimisée en particulier pour Intel Nehalem / Atom , VIA Nanoprocessor , AMD Opteron .
Bibliothèque scientifique GNU: Implémentation multi-plateforme de nombreuses routines numériques. Contient une interface CBLAS.
HP MLIB: Bibliothèque mathématique de HP prenant en charge l'architecture IA-64 , PA-RISC , x86 et Opteron sous HP-UX et Linux .
Intel MKL: La bibliothèque Intel Math Kernel , prenant en charge x86 32 bits et 64 bits, est disponible gratuitement auprès d' Intel . Inclut des optimisations pour les processeurs Intel Pentium , Core et Intel Xeon et Intel Xeon Phi ; prise en charge de Linux , Windows et macOS .
MathKeisan: Bibliothèque mathématique de NEC , prenant en charge l'architecture NEC SX sous SUPER-UX et Itanium sous Linux
Netlib BLAS: L'implémentation de référence officielle sur Netlib , écrite en Fortran 77. [
CBLAS Netlib: Interface de référence C vers le BLAS. Il est également possible (et courant) d'appeler le BLAS Fortran à partir de C.
OpenBLAS: BLAS optimisé basé sur GotoBLAS, prenant en charge les processeurs x86 , x86-64 , MIPS et ARM .
PDLIB/SX: Bibliothèque mathématique du domaine public de NEC pour le système NEC SX-4 .
rocBLAS: Implémentation qui fonctionne sur les GPU AMD via ROCm .
SCSL: La bibliothèque de logiciels de calcul scientifique de SGI contient des implémentations BLAS et LAPACK pour les stations de travail Irix de SGI .
Bibliothèque de performances solaires: BLAS et LAPACK optimisés pour les architectures SPARC , Core et AMD64 sous Solaris 8, 9 et 10 ainsi que Linux.
uBLAS: Une bibliothèque de classes de modèles C++ génériques fournissant des fonctionnalités BLAS. Fait partie de la bibliothèque Boost . Elle fournit des liaisons à de nombreuses bibliothèques accélérées par le matériel dans une notation unifiée. De plus, uBLAS se concentre sur l'exactitude des algorithmes à l'aide de fonctionnalités C++ avancées.

Bibliothèques utilisant BLAS

Tatou: Armadillo est une bibliothèque d'algèbre linéaire C++ qui vise un bon équilibre entre vitesse et facilité d'utilisation. Elle utilise des classes de modèles et possède des liens optionnels vers BLAS/ATLAS et LAPACK. Elle est sponsorisée par NICTA (en Australie) et est sous licence libre.
LAPACK: LAPACK est une bibliothèque d'algèbre linéaire de niveau supérieur basée sur BLAS. Comme BLAS, une implémentation de référence existe, mais de nombreuses alternatives comme libFlame et MKL existent.
Moi: Une bibliothèque numérique générique accélérée par LLVM pour la science et l'apprentissage automatique écrite en D. Elle fournit des sous-programmes d'algèbre linéaire génériques (GLAS). Elle peut être construite sur une implémentation CBLAS.

Bibliothèques similaires (non compatibles avec BLAS)

Élémentaire: Elemental est un logiciel open source pour l'algèbre linéaire directe dense et clairsemée à mémoire distribuée et l'optimisation.
HASEM: est une bibliothèque de modèles C++, capable de résoudre des équations linéaires et de calculer des valeurs propres. Elle est sous licence BSD.
LAMA: La bibliothèque pour les applications mathématiques accélérées (LAMA) est une bibliothèque de modèles C++ pour l'écriture de solveurs numériques ciblant différents types de matériel (par exemple, les GPU via CUDA ou OpenCL ) sur des systèmes de mémoire distribuée , masquant la programmation spécifique au matériel au développeur du programme.
MTL4: La bibliothèque de modèles Matrix version 4 est une bibliothèque de modèles C++ générique offrant des fonctionnalités BLAS clairsemées et denses. MTL4 établit une interface intuitive (similaire à MATLAB ) et une large applicabilité grâce à la programmation générique .

BLAS clairsemé

Plusieurs extensions de BLAS pour la gestion des matrices clairsemées ont été suggérées au cours de l'histoire de la bibliothèque ; un petit ensemble de routines de noyau de matrices clairsemées a finalement été standardisé en 2002.

BLAS par lots

Les fonctions BLAS traditionnelles ont également été portées sur des architectures qui prennent en charge de grandes quantités de parallélisme, comme les GPU . Ici, les fonctions BLAS traditionnelles offrent généralement de bonnes performances pour les grandes matrices. Cependant, lors du calcul, par exemple, de produits matrice-matrice de nombreuses petites matrices à l'aide de la routine GEMM, ces architectures présentent des pertes de performances importantes. Pour résoudre ce problème, en 2017, une version par lots de la fonction BLAS a été spécifiée.

En prenant la routine GEMM ci-dessus comme exemple, la version par lots effectue le calcul suivant simultanément pour plusieurs matrices :

${\displaystyle {\boldsymbol {C}}[k]\leftarrow \alpha {\boldsymbol {A}}[k]{\boldsymbol {B}}[k]+\beta {\boldsymbol {C}}[k]\quad \forall k <semantics> <mrow> <mstyle> <mrow> <mi>C</mi> </mrow> <mo>[</mo> <mi>k</mi> <mo>]</mo> <mo>←</mo> <mi>α</mi> <mrow> <mi>UN</mi> </mrow> <mo>[</mo> <mi>k</mi> <mo>]</mo> <mrow> <mi>B</mi> </mrow> <mo>[</mo> <mi>k</mi> <mo>]</mo> <mo>+</mo> <mi>β</mi> <mrow> <mi>C</mi> </mrow> <mo>[</mo> <mi>k</mi> <mo>]</mo> <mspace></mspace> <mi>∀</mi> <mi>k</mi> </mstyle> </mrow> <annotation>{\displaystyle {\boldsymbol {C}}[k]\leftarrow \alpha {\boldsymbol {A}}[k]{\boldsymbol {B}}[k]+\beta {\boldsymbol {C}}[k]\quad \forall k}</annotation> </semantics> </math><img src=$

L'indice entre crochets indique que l'opération est effectuée pour toutes les matrices d'une pile. Souvent, cette opération est implémentée pour une disposition de mémoire par lots échelonnée où toutes les matrices sont concaténées dans les tableaux , et . ${\displaystyle k <semantics> <mrow> <mstyle> <mi>k</mi> </mstyle> </mrow> <annotation>{\displaystyle k}</annotation> </semantics> </math><img src=$ ${\displaystyle k <semantics> <mrow> <mstyle> <mi>k</mi> </mstyle> </mrow> <annotation>{\displaystyle k}</annotation> </semantics> </math><img src=$ ${\displaystyle A <semantics> <mrow> <mstyle> <mi>UN</mi> </mstyle> </mrow> <annotation>{\displaystyle A}</annotation> </semantics> </math><img src=$ ${\displaystyle B <semantics> <mrow> <mstyle> <mi>B</mi> </mstyle> </mrow> <annotation>{\displaystyle B}</annotation> </semantics> </math><img src=$ ${\displaystyle C <semantics> <mrow> <mstyle> <mi>C</mi> </mstyle> </mrow> <annotation>{\displaystyle C}</annotation> </semantics> </math><img src=$

Les fonctions BLAS par lots peuvent être un outil polyvalent et permettent par exemple une implémentation rapide d' intégrateurs exponentiels et d'intégrateurs Magnus qui gèrent de longues périodes d'intégration avec de nombreux pas de temps. Ici, l' exponentiation matricielle , la partie coûteuse en calcul de l'intégration, peut être implémentée en parallèle pour tous les pas de temps en utilisant des fonctions BLAS par lots.