NumPy-Matrixmultiplikation mit 150 Zeilen C-Code übertroffen

 (salykova.github.io)
1 Punkte von GN⁺ 2024-07-05 | 1 Kommentare | Auf WhatsApp teilen
  • Die Matrixmultiplikation von NumPy stützt sich auf externe BLAS-Bibliotheken, doch diese Implementierung zielt darauf ab, mit reinem C und nur FMA3 sowie AVX2 Single-Thread- und Multi-Thread-Leistung auf BLAS-Niveau zu erreichen
  • Im Zentrum der Performance steht eine Aufteilung von $C$ in kleine Blöcke; ein 16×6-Microkernel wiederholt rank-1 updates in YMM-Registern und reduziert so Speicherzugriffe
  • Bei Matrizen beliebiger Größe wird die Randbehandlung leicht zum Flaschenhals; daher wird eine Kombination aus maskierten Stores und nullaufgefüllerten Puffern verwendet, um den Performanceverlust maskierter Loads zu vermeiden
  • Die Cache-Wiederverwendung wird durch Blocking mit k_c, m_c, n_c sichergestellt; die tatsächliche Spitzenleistung hängt stark vom Tuning von Thread-Anzahl, Kernel-Größe und Tile-Größe ab
  • AVX-512 wurde zugunsten breiterer CPU-Unterstützung ausgeschlossen; auf AVX-512-CPUs kann BLAS daher schneller sein, und auch der Vergleich mit OpenBLAS wurde unter deaktiviertem AVX-512 durchgeführt

Implementierungsziel und Vergleichsgrundlage

  • Der Implementierungscode ist in sgemm.c veröffentlicht und optimiert multithreaded FP32-Matrixmultiplikation auf aktuellen Prozessoren
  • NumPy ist bei linearen Algebra-Operationen wie Matrixmultiplikation auf externe BLAS-Bibliotheken angewiesen
    • Beispiele sind Intel MKL, Accelerate, BLIS, GotoBLAS und OpenBLAS
    • OpenBLAS, GotoBLAS und BLIS sind in C/FORTRAN/Assembly geschrieben und enthalten manuell optimierte Matrixmultiplikations-Implementierungen für verschiedene CPU-Mikroarchitekturen
  • Ziel ist eine Matrixmultiplikations-Implementierung, die ohne Low-Level-Assembly in reinem C geschrieben ist und zugleich die folgenden Bedingungen erfüllt
    • funktioniert für Matrizen beliebiger Größe
    • läuft auf aktuellen x86-64-Prozessoren
    • konkurriert mit bestehenden BLAS-Bibliotheken
    • bleibt einfach und leicht erweiterbar
  • Referenzen sind Simon Boehms Fast Multidimensional Matrix Multiplication on CPU from Scratch, Sergey Slotins Matrix Multiplication, Geohots Can you multiply a matrix? sowie Arbeiten zu GotoBLAS und BLIS

Benchmark-Bedingungen und FLOPS-Berechnung

  • Die Testumgebung besteht aus AMD Ryzen 7 9700X, 32GB DDR5 6000 MHz CL36, OpenBLAS 0.3.26, GCC 13.3 und Ubuntu 24.04.1 LTS
  • Verwendet werden die Compiler-Flags -O3 -march=native -mno-avx512f -fopenmp
  • Für einen fairen Vergleich muss bei der Installation von OpenBLAS ein passendes TARGET gesetzt und AVX-512 deaktiviert werden
    • Zen4/5-Prozessoren werden mit make TARGET=ZEN kompiliert
    • Andernfalls verwendet OpenBLAS standardmäßig AVX-512-Instruktionen
  • Die FP32-Matrixmultiplikation von OpenBLAS wird über die API cblas_sgemm ausgeführt
  • Die Benchmarks verwenden quadratische Matrizen
    • bewertet wird von m=n=k=200 bis m=n=k=10000 in Schritten von 200
    • die Matrixmultiplikation wird n_iter-mal wiederholt, und die mediane Laufzeit dient als Leistungsmaß
  • Werden eine $M \times K$-Matrix $A$ und eine $K \times N$-Matrix $B$ multipliziert, beträgt der gesamte Rechenaufwand $2MNK$ FLOP
    • die Leistung wird als FLOPS=(2*m*n*k)/exec_time berechnet

Theoretische Grenzen und SIMD-Basis

  • Aktuelle x86-64-CPUs verarbeiten mit SIMD-Erweiterungen mehrere Daten parallel
  • Wichtige Instruktionen sind AVX2 und FMA
    • beide verwenden 256-Bit-YMM-Register
    • jedes YMM-Register kann acht 32-Bit-Floats enthalten
  • Die FMA-Instruktion VFMADD231PS führt eine packed-single-Operation in der Form YMM1 = YMM2 * YMM3 + YMM1 aus
  • Auf dem Ryzen 9700X liegt der fused-multiply-add-Durchsatz bei 0.5 cycles/instruction, also zwei Instruktionen pro Takt
  • Theoretisch kann der Ryzen 9700X auf einem einzelnen Kern 32 FLOP pro Takt ausführen
    • die Rechnung lautet 8 floats × 2(add+mul) × 2(1/TP)
    • bei 8 Kernen und einem anhaltenden Takt von 4.7GHz wird der theoretische Multi-Thread-Peak auf 1203 FLOPS geschätzt

Grundlegende Implementierung und Microkernel

  • Die Matrizen werden in column-major-Reihenfolge gespeichert
    • auf einen C-Pointer bezogen wird A[row][col] über ptr[col*M + row] angesprochen
  • Die einfachste Implementierung durchläuft alle Zeilen und Spalten von $C$ und berechnet für jedes Element das Skalarprodukt einer Zeile aus $A$ mit einer Spalte aus $B$
  • Der Kern performanter Implementierungen ist ein Microkernel, der $C$ in Teilmatrizen der Größe $m_R \times n_R$ zerlegt und jede Teilmatrix effizient berechnet
  • Der Kernel initialisiert $\bar{C}$ in Registern mit 0 und iteriert dann entlang der $K$-Dimension
    • lädt den Spaltenvektor von $\bar{A}$ und den Zeilenvektor von $\bar{B}$ in Register
    • berechnet das äußere Produkt der beiden Vektoren und addiert es zu den Akkumulatoren von $\bar{C}$
    • jeder Schritt ist ein rank-1 update
  • Im Vergleich zum Speicherzugriff der naiven Methode von $2K m_R n_R$ reduziert dieser Ansatz die Anzahl der in Register geladenen Elemente auf $(m_R+n_R)K$
  • Da AVX-CPUs 16 YMM-Register besitzen, muss die Kernel-Größe die folgenden Bedingungen erfüllen
    • $(m_R/8) \cdot n_R + m_R/8 + 1 \le 16$
    • $m_R$ muss ein Vielfaches von 8 sein
  • Theoretisch sinken die Speicherzugriffe stärker, je größer und ähnlicher $m_R$ und $n_R$ sind; auf dem Ryzen 9700X liefert jedoch ein 16×6-Kernel die beste Leistung
  • Die Implementierung verwendet Intrinsics aus immintrin.h
    • __m256 ist ein 256-Bit-Vektortyp und repräsentiert den Inhalt eines YMM-Registers
    • _mm256_loadu_ps lädt den Spaltenvektor von A
    • _mm256_broadcast_ss broadcastet einen Skalarwert aus B auf einen Vektor aus acht Floats
    • _mm256_fmadd_ps aktualisiert die Akkumulatoren
    • _mm256_storeu_ps speichert das Ergebnis in den Speicher
  • Der erzeugte Assembler enthält SIMD-FMA-Instruktionen wie vfmadd231ps und vbroadcastss

Padding für Matrizen beliebiger Größe

  • Der grundlegende 16×6-Kernel arbeitet direkt, wenn $M$ und $N$ jeweils Vielfache von 16 beziehungsweise 6 sind
  • Wenn im Randbereich die Spaltenzahl $n$ kleiner als 6 ist, wird die Store-Schleife nur bis j < n ausgeführt
  • Ist die Zeilenzahl $m$ kleiner als 16, speichert _mm256_storeu_ps acht Elemente auf einmal, weshalb maskierte Stores nötig sind
    • _mm256_maskstore_ps speichert nur die Elemente in den Speicher, deren Maskenbits gesetzt sind
    • die Maske wird abhängig von der überlappenden Zeilenzahl $m$ erzeugt
  • Würde auch das Laden im Randbereich mit _mm256_maskload_ps erfolgen, könnte die Kernel-Leistung stark sinken
    • zusätzliche Instruktionen zur Maskenberechnung erzeugen Overhead
    • da $n$ keine Compile-Time-Konstante ist, kann der Compiler die Schleife schwer effizient unrollen
  • Stattdessen wird bei $m \neq m_R$ $\bar{A}$ in einen Puffer kopiert und mit Nullen aufgefüllt; bei $n \neq n_R$ wird auch $\bar{B}$ in einen Puffer kopiert und mit Nullen ergänzt
  • Die zugehörige Implementierung befindet sich in matmul_pad.h

Cache-Blocking und Datenwiederverwendung

  • Zwischen Registern und DRAM liegt eine Hierarchie von CPU-Caches; moderne Desktop-CPUs verwenden typischerweise L1-, L2- und L3-Caches
  • Caches sind schneller als DRAM, aber in der Kapazität begrenzt; daher ist es unmöglich, ganz $A$, $B$ und $C$ gleichzeitig im Cache zu halten
  • Die Matrizen werden in kleine Blöcke zerlegt, die in den Cache geladen werden, und dieselben Daten werden für mehrere rank-1 updates wiederverwendet; dieses Verfahren heißt Cache-Blocking oder Tiling
  • Das Single-Thread-Cache-Blocking hat eine fünffach geschachtelte Schleifenstruktur ähnlich BLIS
    • die äußerste Schleife erzeugt entlang der $N$-Dimension Blöcke $C_j$ und $B_j$
    • die nächste Schleife erzeugt entlang der $K$-Dimension Blöcke $A_j$ und $B_p$
    • $B_p$ wird gepackt zu $\tilde{B}_p$ und bei Bedarf mit Nullen aufgefüllt, um Wiederverwendung im L3-Cache zu erreichen
    • die folgende Schleife erzeugt entlang der $M$-Dimension Blöcke $C_i$ und $A_j$; dabei wird $A_j$ zu $\tilde{A}_j$ gepackt
    • die letzten beiden Schleifen teilen die Cache-Blöcke in Panels der Größen $m_R \times k_c$ und $k_c \times n_R$ und übergeben sie an den Kernel
  • Die gepackten Daten $\tilde{A}_j$ und $\tilde{B}_p$ werden unterschiedlich gespeichert
    • Panels innerhalb von $\tilde{A}_j$ werden in column-major gespeichert
    • Panels innerhalb von $\tilde{B}_p$ werden in row-major gespeichert
  • Die Parameter des Cache-Blockings müssen an die Cache-Größe des jeweiligen CPU-Modells angepasst werden
    • $k_c \times n_c$ ist ein Ausgangspunkt zum Füllen des L3-Caches
    • $m_c \times k_c$ ist ein Ausgangspunkt zum Füllen des L2-Caches
    • $k_c \times n_R$ ist ein Ausgangspunkt zum Füllen des L1-Caches
  • In der Praxis liefern oft größere Werte als die theoretischen besseren Durchsatz, und da die CPU die Cache-Platzierung selbst verwaltet, müssen Schleifen und Zugriffsmuster auf Algorithmus-Ebene entworfen werden
  • Die Implementierung befindet sich in matmul_cache.h

Feintuning des Kernels

  • Statt Akkumulatoren als Array wie __m256 C_buffer[6][2] zu definieren, werden die Akkumulatorvariablen explizit ausgerollt deklariert
  • Dieser Ansatz hilft GCC, den Code besser zu optimieren und Register-Spilling zu vermeiden
  • Auch die Maskenberechnung wurde auf Vektor-Instruktionen umgestellt
    • dazu wird ein statisches Array mask[32] verwendet sowie _mm256_cvtepi8_epi32 und _mm_loadu_si64
  • Die betreffende Implementierung steht in matmul_micro.h

Multithreading-Strategie

  • Parallelisiert werden sowohl arithmetische Operationen als auch das Packing
  • Die 5., 4. und 3. Schleife außerhalb des Microkernels iterieren in Einheiten der Cache-Blockgröße
    • um alle Threads auszulasten, muss die Anzahl der Iterationen mindestens so groß wie die Thread-Anzahl sein
    • die Dimensionen der Eingabematrix sollten ungefähr mindestens Thread-Anzahl × Cache-Blockgröße betragen
  • Auf einem Ryzen 9700X mit einem einzelnen Thread zeigten Cache-Blockgrößen von $n_c=1535$, $m_c=1024$ gute Leistung
    • um alle 8 Kerne auszulasten, wäre mindestens eine Dimension von $\max(m_c,n_c) \times 8 = 1535 \times 8 = 12280$ nötig
  • Umgekehrt wiederholen die letzten beiden Schleifen kleine Blöcke von $m_R$ und $n_R$ und eignen sich daher gut für Parallelisierung
    • typischerweise sind $m_R$ und $n_R$ kleiner als 20
    • wählt man $m_c$ und $n_c$ als Vielfache der Kernanzahl, lässt sich die Arbeit gleichmäßig verteilen
  • Auf dem Ryzen 9700X liefert das gemeinsame Parallelisieren der beiden inneren Schleifen mit #pragma omp parallel for collapse(2) num_threads(NTHREADS) die beste Leistung
  • Auf Prozessoren mit vielen Kernen, besonders mit mehr als 16 Kernen, kann man verschachtelte Parallelität und die Parallelisierung von 2 bis 3 Schleifen in Betracht ziehen
  • Auch das Packing von $\tilde{A}$ und $\tilde{B}$ wird mit OpenMP parallelisiert
    • pack_blockA wird parallelisiert, indem über mc in Schritten von MR iteriert wird
    • pack_blockB wird parallelisiert, indem über nc in Schritten von NR iteriert wird
  • Für die Multi-Thread-Implementierung zeigten auf dem Ryzen 9700X die folgenden Parameter gute Leistung
    • $m_c = m_R \times \text{Anzahl der Threads} \times 5$
    • $n_c = n_R \times \text{Anzahl der Threads} \times 50$
  • Die finale Multi-Thread-Implementierung befindet sich in matmul_parallel.h

Verwandte Beiträge

1 Kommentare

 
GN⁺ 2024-07-05
Meinungen auf Hacker News

  • Wenn die Kernaussage dieses Beitrags ist, dass im Allgemeinen noch Performance-Potenzial übrig ist, dann unterschätzt er das Ausmaß der möglichen Verbesserungen eher. Und das, obwohl in Bibliotheken für Matrixmultiplikation deutlich mehr Aufwand steckt als in der meisten Software.
    Wenn Code nicht bereits stark optimiert ist, sind Verbesserungen um den Faktor 10 bis über 1000 gegenüber bestehendem Code oft ohne großen Aufwand möglich. Grob nach Wichtigkeit sortiert: Am wichtigsten ist, ob der gewählte Algorithmus passend ist und ob sich die Arbeit selbst vermeiden lässt; ebenfalls viel bringt es, schwere Operationen wie Kernel-Roundtrips oder malloc zu reduzieren.
    Bei der Vektorisierung sind explizite Vektor-Intrinsics gut, aber oft reicht es schon, Daten von einem Array von Strukturen in Strukturen/Arrays von Arrays umzubauen, damit derselbe Maschinencode entsteht. Cache-Effizienz ist ebenfalls wichtig, und bei parallelem Code wird es komplizierter, wenn Daten nicht sauber pro Thread isoliert sind, etwa bei False Sharing. Schließlich sind auch hardware-spezifische Optimierungen möglich, etwa mit Intrinsics oder handgeschriebenem Assembly.

    • Auch den Einfluss des Netzwerks darf man nicht vergessen. Ich habe einmal eine verteilte Abfrage gefunden, die rund eine Million Zeilen über das Netzwerk holte und sie dann jointe, sodass nur 5 bis 10 Zeilen übrig blieben; durch die Änderung ließ sich die Performance um mehrere hundertfach verbessern.
      Nachdem die Abfrage so geändert wurde, dass der Join auf dem entfernten Server stattfand und nur noch 5 bis 10 Zeilen über das Netzwerk gesendet wurden, war sie sofort schneller. Fixe Overheads und Latenz gibt es zwar immer, aber wenn man weit mehr Daten als nötig über eine Netzwerkverbindung schickt, ruiniert das am Ende die Performance. Lesenswert zum Einfluss von Latenz ist auch „It's the latency, stupid“: http://www.stuartcheshire.org/rants/latency.html
      Insgesamt stimme ich den obigen Punkten und ihrer groben Reihenfolge zu.
    • „Ist die Algorithmuswahl passend?“ ist in der Praxis teilweise zu Cargo Cult geworden. Ein „schnellerer“ Algorithmus kann schreckliche konstante Faktoren haben, sodass die Variante, die mehr Arbeit erledigt, oft tatsächlich besser performt.
      Viele Vorstellungsgespräche sind statt einer Prüfung, ob man herleiten, benchmarken und beheben kann, warum eine Implementierung langsam ist, zu Quizzen über obskure Algorithmen zum Auswendiglernen geworden – nach dem Motto „weil Google das so macht“.

  • Gängige Coding-Patterns sind nicht ausreichend auf die Hardware spezialisiert und lassen daher viel Performance liegen. Dieser Beitrag ist ein interessantes Beispiel dafür; eine weitere klassische Demonstration ist „There's plenty of room at the top“.
    https://www.science.org/doi/10.1126/science.aam9744

  • Um das zu verstehen, sind die Papers im BLIS-Repository praktisch die Referenz. Ich weiß nicht, warum man glauben sollte, optimiertes BLAS liefere keine Performance; bei ausreichend großen Matrizen sollte man mehr als 90 % der CPU-Spitzenleistung erwarten.
    Als ich zuletzt nachgesehen habe, lag serielles OpenBLAS im Großen und Ganzen auf MKL-Niveau, und BLAS implementiert als grundlegenden Baustein der linearen Algebra nicht matmul, sondern GEMM. Ich verstehe auch nicht, warum man numpy statt eines Benchmark-Frameworks verwendet; auf Zen sollte man meiner Meinung nach mit AMDs BLAS vergleichen, also mit der BLIS-basierten Implementierung. BLIS hatte früher die bessere Geschichte bei der Parallelisierung als OpenBLAS, und AMD BLIS hat auch eine Umschaltung auf Implementierungen für „kleine“ Dimensionen; ob es das derzeit in OpenBLAS gibt, weiß ich nicht.
    Für die Vektorisierung von Mikro-Kernels braucht man nicht zwingend SIMD-Intrinsics; ein guter C-Compiler vektorisiert vollständig und rollt die Schleifen auch aus. Die reinen C-Mikro-Kernels von BLIS erreichen bei passenden Blockgrößen mehr als 80 % der Performance der handoptimierten Haswell-Implementierung. Der Unterschied liegt vermutlich am Prefetching, aber ganz genau verstehe ich das nicht.

    • SIMD-Intrinsics und manuelles Loop-Unrolling sind ganz klar nötig. Genau deshalb vektorisieren und entrollen alle BLAS-Bibliotheken ihre Schleifen von Hand.
      Auch moderne Compiler schaffen Auto-Vektorisierung und Loop-Unrolling nicht zuverlässig und korrekt mit 100 % Erfolgsquote.

  • Der Artikel und die Implementierung sehen gut aus, aber ich frage mich, worin das „Geheimnis“ liegt. OpenBLAS wurde über Jahrzehnte hinweg genau für dieses Problem in Assembly+C optimiert – wie kann man das schlagen?
    Es geht ausführlich um Caching usw.; ich frage mich, ob BLAS solche Dinge nicht nutzt oder ob es einfach besser auf einen bestimmten Prozessor zugeschnitten ist.

    • OpenBLAS ist nicht unbedingt so stark auf bestimmte moderne Architekturen optimiert. Die Matrizen waren auch nicht besonders groß, und bei numpy gibt es cffi-Overhead.
      Der Leistungsunterschied war beim Peak-Durchsatz viel ausgeprägter als beim durchschnittlichen Durchsatz, aber es gibt kaum Anwendungen, bei denen der Peak entscheidend ist. Der gezeigte Benchmark-Code scheint auf der numpy-Seite durch den Python-Allocator zu gehen, während die C-Implementierung den Allocator nicht durchläuft; das wäre also die erste Stelle, an der man nach Fehlern oder Inkonsistenzen im Mikrobenchmark suchen sollte. Viele numpy-Routinen unterstützen In-place-Operationen, daher sollte man für beide Seiten ausdrücklich die In-place-Versionen benchmarken.
      numpy hat außerdem unabhängig von der zugrunde liegenden Implementierung Bounds-Checks und Fehlerbehandlung; das ist ein Grund, warum es bei kleinen Matrizen sogar deutlich langsamer sein kann als reine Python-Listen. Wenn ein paar Tausend Zyklen reiner Overhead dazukommen, wird es schwer, schnell zu sein.
      Diese Implementierung ist ein ziemlich prinzipientreuer Ansatz, die relevanten Caches auszulasten, und in gewisser Weise naheliegend; aber klare Engineering-Verbesserungen verdienen es, in solchen Diskussionen hervorgehoben zu werden. Auch in OpenBLAS ist viel Arbeitskraft geflossen, aber es ist unwahrscheinlich, dass dort an alles gedacht wurde. Um das sauber zu erklären, bräuchte es eine tiefgehende Analyse beider Codebasen.
    • OpenBLAS zu schlagen ist weder überraschend noch beispiellos. Zum Beispiel hat das die Lineare-Algebra-Bibliothek Mir der Sprache D schon vor einigen Jahren geschafft [1].
      Für C++- und C-Implementierungen siehe die Metaprogramming-Ansätze [2], [3]. Wirklich überraschend ist, dass viele moderne Sprachen wie Matlab, Julia und Mojo immer noch auf OpenBLAS setzen – natürlich werden sie jeweils ihre Gründe haben.
      [1] Numeric age for D: Mir GLAS is faster than OpenBLAS and Eigen (2016):
      http://blog.mir.dlang.io/glas/benchmark/openblas/2016/09/23/...
      [2] Vastly outperforming LAPACK with C++ metaprogramming (2018):
      https://wordsandbuttons.online/vastly_outperforming_lapack_w...
      [3] Outperforming LAPACK with C metaprogramming (2018):
      https://wordsandbuttons.online/outperforming_lapack_with_c_m...
    • -march=native kann einen Vorteil bringen, weil damit für das genaue CPU-Modell kompiliert wird. numpy wurde sehr wahrscheinlich für ein allgemeineres und älteres x86-64-Ziel kompiliert.
      Auf einer Ryzen-CPU nutzt -march=native vermutlich v4, während numpy eher v1 oder v2 als Ziel haben dürfte.
      https://en.wikipedia.org/wiki/X86-64#Microarchitecture_level...
    • numpy 2.0 integriert Google Highway, um SIMD auf mehreren Mikroarchitekturen besser zu nutzen; der Vergleich dürfte für numpy dadurch besser ausfallen.

  • Der Artikel ist gut, und es ist großartig, dass die Benchmarks leicht reproduzierbar gemacht wurden. Auf meinem 16-Core Xeon W-2245 3.90GHz erledigte matmul.c eine 8192x8192-Matrixmultiplikation mit gcc -O3 in 1,41 Sekunden und mit clang -O2 in 1,47 Sekunden; NumPy brauchte 1,07 Sekunden.
    Mit einem AVX-512-Kernel dürfte es deutlich schneller werden. Ein weiterer Grund für die enttäuschende Performance könnte OpenMP sein; meiner Erfahrung nach kann man Overhead reduzieren, wenn man einen Thread-Pool mit pthreads explizit verwaltet. Statt die CPU-Anzahl hart zu codieren, wäre außerdem sysconf(_SC_NPROCESSORS_ONLN) besser.

  • Es gibt keinen Grund, die Last ungleich zu verteilen, indem die eine Seite Python und die andere C ist. Man hätte einen Apples-to-Apples-Vergleich machen können, bei dem beide Seiten in C geschrieben sind: die eine ruft eine BLAS-Bibliothek auf, die andere diese Implementierung.

    • Hier ist der Vergleich mit Python richtig. Denn die heute populärste Art, solche Berechnungen durchzuführen, ist Python mit numpy.
      Der Overhead ist nicht riesig, aber wie auch anderswo in diesem Thread gesagt wurde, ist es wichtig, korrekt aufzurufen. Naiven numpy-Code gegen abgestimmten C-Code antreten zu lassen, ist eindeutig kein fairer Vergleich.

  • Auch wenn es kein Hot Path ist: Die Ineffizienz bei der Maskenerzeugung stört mich, nämlich die Verwendung von bit_mask. Effizientere Varianten wären, ein globales Konstanten-Array der Form {-1,-1,...,0,0,...} anzulegen und ab den Element-Offsets 16-m bzw. 8-m zu laden, oder einen konstanten Vektor {0,1,2,3,4,...} mit dem gebroadcasteten m und m-8 zu vergleichen.
    Allerdings betrifft das nur eine Spalte der Matrix, und die anschließende maskload/maskstore-Schleife dauert deutlich länger, daher ist das nur eine sehr kleine Nörgelei. Vor allem Stores sind selbst auf Zen 4 immer noch langsam[1], und AVX-512-Instruktionen sind 6-mal schneller, obwohl der einzige Unterschied darin besteht, dass sie die Maske aus Maskenregistern beziehen. Clang vektorisiert die Shifts ohnehin automatisch, es dürfte also nur etwa 2- bis 3-mal langsamer sein als mein Vorschlag.
    [1]: https://uops.info/table.html?search=vmaskmovps&cb_lat=on&cb_...

    • Ich bin der Autor. Optimierung von C-Code und die Verwendung von Intrinsics sind für mich wirklich Neuland, daher bin ich kein Experte auf diesem Gebiet, möchte aber mehr lernen.
      Ich weiß Feedback, das neue Perspektiven eröffnet, wirklich zu schätzen. „Ein konstantes globales Array anlegen und daraus laden“ war meiner Erinnerung nach in meinen Tests etwas langsamer als das Verschieben der Bitmaske, aber ich werde es noch einmal testen, um sicherzugehen. Die Variante, „einen konstanten Vektor {0, 1, 2, 3, 4, ...} mit dem gebroadcasteten m und m-8 zu vergleichen“, ist eine gute Idee; ich werde sie ausprobieren.
    • Wenn man ein globales Konstanten-Array anlegt, kann man die Elemente als int8_t ablegen und beim Laden von Byte auf int32_t vorzeichenerweitern. Die Kombination _mm_loadu_si64 / _mm256_cvtepi8_epi32 dürfte zu einer einzelnen vpmovsxbd-Instruktion mit Speicheroperand kompiliert werden.
      Damit passt das gesamte Konstanten-Array in eine Cache Line, wenn es mit alignas(32) korrekt ausgerichtet ist. Im Anwendungsfall des Originals werden zwei Masken benötigt, daher ist die zweite vpmovsxbd-Instruktion mit Sicherheit ein L1D-Cache-Hit, was gut passt.

  • Wie wäre es mit jarts tinyBLAS?
    https://hacks.mozilla.org/2024/04/llamafiles-progress-four-m...
    Und https://justine.lol/matmul/

    • Ich habe gestern ausführlich mit Justine gesprochen, und auf dieser Workstation scheint diese Implementierung mindestens 2-mal schneller zu sein als tinyBLAS. Die vollständige Diskussion ist im Mozilla-AI-Discord: https://discord.com/invite/NSnjHmT5xY

  • Abgesehen von Benchmarks: Warum sollte man die Matrixmultiplikation selbst multithreaden? Wäre es in der Praxis nicht vorteilhafter, Multithreading auf der Ebene des Algorithmus einzusetzen, der die Multiplikation verwendet?

    • In HPC macht man das tatsächlich meist so. Allerdings kann es bestimmten Arten von R-Code schon leicht helfen, einfach auf ein paralleles BLAS umzusteigen.
      HPC-Code hat jedoch meist GEMM nicht als Bottleneck.

  • Ich habe bisher nur darübergeschaut, aber der Artikel enthält viele Details und Erklärungen. Er wirkt wie ein ziemlich guter Text darüber, wie schnelle Matrixmultiplikation unter Berücksichtigung von Architekturüberlegungen implementiert wird; ich habe ihn auf meine Leseliste gesetzt.