Computational Sciences Center

Vektorisierung

Viele heute verfügbare Prozessoren bieten in begrenztem Umfang die Möglichkeit, Operationen auf Vektoren aus mehreren Werten statt auf einzelne Werte anzuwenden. Für numerische Programme sind dabei vor allem die Befehlssätze SSE und AVX von Interesse, die mit 128- oder 256bittigen Operanden arbeiten.

Eine knappe Einführung in SSE findet sich in Kapitel 4 meines HPC-Skripts. Für eine Übersicht über die verfügbaren Intrinsics sei auf eine Datenbank der Firma Intel verwiesen.

SSE

Moderne x86-Prozessoren bieten in der Regel Unterstützung für die Streaming SIMD Extension (SSE). Sie besteht aus einer Reihe von Maschinenbefehlen, die mit 128bittigen Vektorregistern arbeiten. Falls der verwendete Compiler nicht dazu in der Lage ist, den vorliegenden Quelltext automatisch zu vektorisieren, können wir ihm unter die Arme greifen, indem wir die in "xmmintrin.h" (SSE) beziehungsweise "emmintrin.h" (SSE2) definierten Intrinsics verwenden. Dabei handelt es sich um Funktionen, die der Compiler automatisch in korrespondierende Maschinenbefehle umsetzt, so dass wir im Wesentlichen auf Maschinenebene programmieren.

Die folgenden Dateien lassen sich mit einem Makefile automatisch übersetzen. Für die Zeitmessung sind die Dateien stopwatch.c und stopwatch.h erforderlich.

Arithmetische Operationen: axpy

Die aus BLAS bekannte axpy-Operation berechnet die Linearkombination zweier Vektoren, addiert also zu einem gegebenen Vektor ein Vielfaches eines anderen Vektors hinzu.

Mit SSE können wir diese Operation realisieren, indem wir die Intrinsics _mm_mul_ps für die komponentenweise Multiplikation und _mm_add_ps für die komponentenweise Addition verwenden. Dabei steht "ps" für "packed single precision", ein 128bittiges Vektorregister wird als Vektor aus vier 32bittigen IEEE-Gleitkommazahlen einfacher Genauigkeit interpretiert.

Solche Daten können mit _mm_load_ps und _mm_store_ps aus dem Speicher gelesen und in den Speicher geschrieben werden. Ein Sonderfall ist _mm_load1_ps, bei dem eine einzige Zahl in alle Komponenten des Vektorregisters geschrieben wird. Wir verwenden diese Variante, um den Skalierungsparameter zu vervielfältigen.

sse_axpy.c

Die GNU Compiler Collection bietet eine sehr elegante Alternative zu den Intrinsics: Sie unterstützt unmittelbar Vektorregister und definiert, wie bestimmte arithmetische Operationen auf diese Register ausgedehnt werden. Für einfache Aufgaben wie die vorliegende können wir damit unser Programm erheblich abkürzen, können es dann aber nur noch mit dem GNU-Compiler übersetzen.

gcc_axpy.c

Gleitkomma- und ganze Zahlen: exp

Gerade bei langen Vektoren werden wir keine sehr beeindruckende Beschleunigung bei der axpy-Operation beobachten, weil die verfügbare Speicherbandbreite den Geschwindigkeitsgewinn begrenzt. Deshalb untersuchen wir eine Operation, die einen erheblich höheren Rechenaufwand pro Operand aufweist, nämlich die Berechnung der Exponentialfunktion.

Da SSE dafür keinen Maschinenbefehl zur Verfügung stellt, basteln wir uns selbst eine Approximation der Exponentialfunktion. Als Ausgangspunkt dient uns die konventionelle Exponentialreihe, die wir mit Hilfe des Horner-Schemas auswerten. Allerdings konvergiert diese Reihe für große Operanden sehr langsam, so dass wir einen Trick anwenden: Wir zerlegen den Operanden in ein ganzzahliges Vielfaches von ln(2) und einen Rest, der zwischen 0 und ln(2) liegt. Für diesen Rest konvergiert die Exponentialreihe hinreichend schnell. Für das ganzzahlige Vielfach haben wir exp(k ln(2)) = 2^k, so dass wir durch Manipulation des Exponenten der IEEE-Gleitkommadarstellung das Ergebnis der Exponentialreihe mit dieser Zweierpotenz multiplizieren können. In C99 gibt es dafür die Funktion ldexp, für SSE müssen wir explizit auf die Binärdarstellung der Gleitkommazahl zurückgreifen. Dabei helfen uns die in SSE2 definierten Intrinsics _mm_cvttps_epi32 (Abrunden zu einer 32bittigen Ganzzahl), _mm_vvtepi32_ps (Konvertieren einer 32bittigen Ganzzahl in eine 32bittige Gleitkommazahl), _mm_slli_epi32 (Bit-Shift nach links), _mm_castsi128_ps (Interpretation eines Vektorregisters als 128bittige Ganzzahl) und _mm_castpi_si128 (Interpretation eines Vektorregisters als vier 32bittige Gleitkommazahlen). Die letzten beiden Intrinsics führen keine Operationen aus, sie signalisieren dem Compiler nur, dass wir ein Vektorregister unterschiedlich interpretieren wollen.

sse_exp.c

Datenabhängige Fallunterscheidungen: if

Vektorrechner haben traditionell Schwierigkeiten mit Fallunterscheidungen, die von den Daten eines Vektorregisters abhängen: Ein Vektorrechner führt grundsätzlich immer eine Operation für alle Komponenten eines Vektorregisters aus, also erfordert es etwas Arbeit, um unterschiedliche Operationen für unterschiedliche Komponenten zustande zu bringen. SSE und AVX verfolgen dabei einen einfachen Ansatz: Mit Funktionen wie _mm_cmplt_ps wird eine Bitmaske erzeugt, in der diejenigen Komponenten gesetzt sind, in denen eine Bedingung zutrifft, in unserem Fall ist das erste Argument kleiner (less than, deshalb cmplt) als das zweite ist. Bei einer Fallunterscheidung führen wir einfach beide Zweige aus und verwenden anschließend die Bitmaske und bitweise logische Operationen wie _mm_and_ps, _mm_andnot_ps und _mm_or_ps, um aus den Ergebnissen beider Zweige das Endergebnis zusammenzustellen.

Im Beispiel berechnen wir das Maximum. Dabei gibt es zwei unterschiedliche Zugänge: Wir können ausschließlich mit logischen Operationen arbeiten, wir können aber auch eine Mischung aus logischen und arithmetischen Operationen verwenden. Da die bitweise 0 auch die Gleitkomma-Null ist, können wir das Maximum von x und y berechnen, indem wir y-x zu x addieren, falls x kleiner als y ist, und sonst null addieren. Die logisch-arithmetische Variante ist auf manchen Prozessoren etwas schneller. Das Beispiel ist allerdings etwas ungeschickt gewählt, denn es gibt für die Berechnung des Maximums auch das Intrinsic _mm_max_ps.

sse_if.c

Umsortieren von Komponenten: addmul

Gelegentlich wollen wir die Komponenten eines Vektorregisters umsortieren. Dafür können wir das Intrinsic _mm_shuffle_ps einsetzen, das aus den Komponenten zweier Vektorregistern ein neues Vektorregister zusammenstellt. Dabei legt eine 8bittige Zahl fest, welche Komponenten wohin kopiert wird. Diese Zahl wird in vier 2bittige Zahlen zerlegt, die jeweils angeben, welche Komponente übernommen werden soll.

Als Beispiel betrachten wir die Berechnung der Multiplikation eines vierdimensionalen Vektors mit einer vierdimensionalen quadratischen Matrix. Die Operation wird spaltenweise durchgeführt: Zuerst wird die erste Spalte mit der ersten Komponente des Eingabevektors multipliziert und zu dem Ergebnis addiert, dann die zweite, dritte und vierte. _mm_shuffle_ps kommt zum Einsatz, um jeweils eine Komponente des Eingabevektors auszuwählen.

sse_addmul.c

AVX

Viele moderne x86-Prozessoren bieten neben SSE mit seinen 128bittigen Vektorregistern auch Unterstützung für AVX (Advanced Vector Extension) mit 256bittigen Registern. Ein derartiges Register kann acht Gleitkommazahlen einfacher oder vier doppelter Genauigkeit aufnehmen, und Rechenoperationen mit diesen Registern erfordern in vielen Fällen dieselbe Zeit wie ihre 128bittigen Gegenstücke. Allerdings gibt es Ausnahmen: Manche Prozessoren reduzieren die Taktfrequenz, sobald AVX-Befehle eingesetzt werden. Andere bieten keine echte 256bittigen Rechenwerke, sondern simulieren sie lediglich. In diesen Fällen sind AVX-Programme nicht so schnell, wie man es auf den ersten Blick erwarten würde. Auch AVX-Befehle können über Intrinsics verwendet werden, die in der Headerdatei "immintrin.h" definiert sind.

Die folgenden Beispielprogramme können mit diesem Makefile übersetzt werden und greifen ebenfalls auf stopwatch.c und stopwatch.h für die Zeitmessung zurück.

Arithmetische Operationen: axpy

Die klassische axpy-Operationen sieht ihrem SSE-Gegenstück sehr ähnlich, wir setzen lediglich _mm256_add_ps und _mm256_mul_ps für Additionen und Multiplikationen ein und verwenden _mm256_broadcast_ss, um den Skalierungsparameter in alle Komponenten eines Vektorregisters zu übernehmen. Natürlich müssen wir auch darauf achten, die Schrittweite für die for-Schleife anzupassen.

avx_axpy.c

Gleitkomma- und ganze Zahlen: exp

Bei der Berechnung der Exponentialfunktion können wir ebenfalls wie im SSE-Fall vorgehen, allerdings ergibt sich eine kleine Schwierigkeit, weil AVX nicht alle Ganzzahloperationen bietet, die wir brauchen, um den Exponenten der Gleitkommazahl zu manipulieren. Deshalb sind wir dazu gezwungen, mit _mm256_extractf128_ps und _mm256_insertf128_ps die 256bittigen Vektorregister in zwei 128bittige zu zerlegen, auf die wir die SSE2-Ganzzahloperationen anwenden können. Mit AVX2 könnten wir das Programm erheblich vereinfachen, denn diese Erweiterung bietet die nötige Operation _mm256_slli_epi32.

avx_exp.c

Fallunterscheidungen und Wurzel: sincos

Fallunterscheidungen lassen sich ebenfalls wie mit SSE umsetzen: _mm256_cmp_ps berechnet eine Bitmaske, die das Ergebnis einer Vergleichsoperation darstellt. Anders als bei SSE wird der durchzuführende Vergleich dabei über einen separaten Parameter festgelegt. Anschließend können wir mit _mm256_and_ps, _mm256_andno_ps und _mm256_or_ps aus den für die verschiedenen Fälle berechnen Teilergebnissen das Gesamtergebnis zusammensetzen.

Bei der Berechnung des Sinus kommen allerdings noch einige Feinheiten hinzu: Da die verwendete Taylor-Entwicklung nur in der Nähe der Null schnell konvergiert, verwenden wir _mm256_floor_ps, um ein Vielfaches der Periode zu subtrahieren und so einen Operanden zwischen -pi und pi zu erhalten. Anschließend spiegeln wir an -pi/2 oder pi/2, um das Argument auf das Intervall von -pi/2 bis pi/2 einzuschränken. Nun können wir eine konventionelle Taylor-Entwicklung einsetzen, bei der wir natürlich ausnutzen, dass nur Terme zu ungeraden Ableitungen auftreten.

Wir können mit dem Satz von Pythagoras auch einfach aus dem Sinus den Cosinus rekonstruieren. Dazu müssen wir uns lediglich an geeigneter Stelle die entsprechenden Vorzeichen merken und können den Algorithmus sogar noch etwas verfeinern: Auf dem Intervall von -pi/4 bis pi/4 approximieren wir den Sinus, ansonsten approximieren wir den Cosinus. Für die Rekonstruktion der jeweils anderen Winkelfunktion verwenden wir dann Pythagoras.

Dafür brauchen wir die Wurzelfunktion. Die eingebaute Fassung _mm256_sqrt_ps ist allerdings unter Umständen nicht besonders effizient, so dass wir an ihrer Stelle lieber die reziproke Wurzelfunktion 1/sqrt(x) approximieren, aus der sich dann einfach die Wurzel gewinnen lässt. Für die reziproke Wurzel können wir mit _mm256_rsqrt_ps sehr schnell eine Näherung erhalten, allerdings eine schlechte. Ein Schritt des Newton-Verfahrens bringt die gewünschte Genauigkeit und lässt sich glücklicherweise nur mit Multiplikationen und einer Subtraktion umsetzen, also mit Grundrechenarten, die der Prozessor sehr schnell ausführen kann.

avx_sincos.c