Ist es möglich, vectorize myNum + = a [b [i]] * c [i]; auf x86_64?

Question

Was intrinsics würde ich vektorisieren die folgende (wenn es überhaupt möglich ist zu vektorisiert) auf der x86_64?

double myNum = 0;
for(int i=0;i<n;i++){
    myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double
}

Answer 1

Hier ist mein gehen an ihm, vollständig optimiert und getestet:

#include <emmintrin.h>

__m128d sum = _mm_setzero_pd();
for(int i=0; i<n; i+=2) {
    sum = _mm_add_pd(sum, _mm_mul_pd(
        _mm_loadu_pd(c + i),
        _mm_setr_pd(a[b[i]], a[b[i+1]])
    ));
}

if(n & 1) {
    sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1]));
}

double finalSum = _mm_cvtsd_f64(_mm_add_pd(
    sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1))
));

Das erzeugt sehr schönen Assembler-Code mit gcc -O2 -msse2 (4.4.1).

Wie Sie sagen, eine noch n hat, wird diese Schleife schneller gehen sowie eine ausgerichtete c. Wenn Sie für eine noch schnellere Zeiten Ausführung c, ändern _mm_loadu_pd zu _mm_load_pd auszurichten.

Answer 2

würde ich durch Abrollen der Schleife starten. So etwas wie

double myNum1 = 0, myNum2=0;
for(int i=0;i<n;i+=2)
{
    myNum1 += a[b[ i ]] * c[ i ];
    myNum2 += a[b[i+1]] * c[i+1];
}
// ...extra code to handle the remainder when n isn't a multiple of 2...
double myNum = myNum1 + myNum2;

Hoffentlich, die die Compiler die Lasten mit dem arithmetischen verschachteln ermöglicht; Profil und Blick auf die Montage zu sehen, ob eine Verbesserung gibt. Im Idealfall wird die Compiler-SSE-Befehle generieren, aber ich bin nicht, ob das in der Praxis geschieht.

Abrollen weiter könnte lassen Sie dies tun:

__m128d sum0, sum1;
// ...initialize to zero...
for(int i=0;i<n;i+=4)
{
    double temp0 = a[b[ i ]] * c[ i ];
    double temp1 = a[b[i+1]] * c[i+1];
    double temp2 = a[b[i+2]] * c[i+2];
    double temp3 = a[b[i+3]] * c[i+3];
    __m128d pair0 = _mm_set_pd(temp0, temp1);
    __m128d pair1 = _mm_set_pd(temp2, temp3);
    sum0 = _mm_add_pd(sum0, pair0);
    sum1 = _mm_add_pd(sum1, pair1);
}
// ...extra code to handle the remainder when n isn't a multiple of 4...
// ...add sum0 and sum1, then add the result's components...

(Entschuldigung für die Pseudo-Code am Anfang und Ende, ich meine, der wichtige Teil der Schleife war). Ich weiß nicht sicher, ob dies schneller sein wird; es hängt von den verschiedenen Latenzen und wie gut die Compiler alles neu anordnen können. Stellen Sie sicher, dass Sie das Profil vor und nach, um zu sehen, ob es eine tatsächliche Verbesserung.

Ich hoffe, das hilft.

Answer 3

Intel-Prozessoren können zwei Gleitkomma-Operationen ausgeben, sondern eine Last pro Zyklus, so Speicherzugriffe die engste Beschränkung ist. Mit diesem Hintergrund sollte ich zuerst verpackt Lasten verwenden, um die Anzahl der Ladeanweisungen zu reduzieren und verwendete Arithmetik verpackt, nur weil es praktisch war. Ich habe, dass die Speicherbandbreite zu sättigen, da realisierte das größte Problem sein, und alle mit SSE-Befehlen des Herumspiel könnte vorzeitige Optimierung, wenn der Punkt gewesen war gehen Sie den Code ein, um schnell und nicht lernen vektorisieren.

SSE

Die wenigsten möglichen Belastungen ohne Annahme auf den Indizes in b erfordert Abrollen der Schleife viermal. Eine 128-Bit-Last erhält vier Indizes von b, zwei 128-Bit-Lasten erhalten jeweils ein Paar benachbarter verdoppelt mich von c, und dem Sammeln von a erforderlich unabhängigen 64-Bit-Lasten. Das ist ein Stockwerk von 7 Zyklen pro vier Iterationen für Seriencode. (Genug, um meine Speicherbandbreite zu sättigen, wenn der Zugang zu a nicht schön zwischenspeichert). Ich habe ein paar ärgerliche Dinge ausgelassen wie eine Anzahl von Iterationen Handhabung, die nicht ein Vielfaches von 4 ist.

entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c)
  xorpd xmm0, xmm0
  xor r8, r8
loop:
  movdqa xmm1, [rdx+4*r8]
  movapd xmm2, [rcx+8*r8]
  movapd xmm3, [rcx+8*r8+8]
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm4, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm4, [rsi+8*r10]
  punpckhqdq xmm1, xmm1
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm5, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm5, [rsi+8*r10]
  add    r8,   4
  cmp    r8,   rdi
  mulpd  xmm2, xmm4
  mulpd  xmm3, xmm5
  addpd  xmm0, xmm2
  addpd  xmm0, xmm3
  jl loop

Immer die Indizes aus ist der komplizierteste Teil. movdqa Lasten 128 Bits von ganzzahligen Daten von einer 16-Byte-ausgerichteten Adresse (Nehalem Latenz Strafen hat zum Mischen der „integer“ und „float“ SSE-Befehle). punpckhqdq bewegt hohe 64 Bits auf 64 Bits niedriger, aber in der Integer-Modus im Gegensatz zu der mehr einfach benannt movhlpd. 32-Bit-Verschiebungen sind in den Mehrzweckregistern erfolgen. movhpd lädt ein Doppel in den oberen Teil eines XMM-Register ohne den unteren Teil zu stören -. diese verwendet werden, um die Elemente von a in gepackte Register direkt zu laden

Dieser Code deutlich schneller als der obige Code, der als der einfache Code wiederum schneller ist, und auf jedem Zugriffsmuster, aber der einfache Fall B[i] = i wo die naive Schleife tatsächlich schnellste ist. Ich habe auch versucht, ein paar Sache wie eine Funktion um SUM(A(B(:)),C(:)) in Fortran, die im wesentlichen äquivalent zu der einfachen Schleife endeten.

I getestet auf einem Q6600 (65 nm Core 2 bei 2,4 GHz) mit 4 GB DDR2-667 Speicher, in 4 Module. Testen Speicherbandbreite gibt etwa 5333 MB / s, so scheint es, als ob ich nur einen einzigen Kanal zu sehen. Ich bin Kompilieren mit Debians gcc 4.3.2-1.1, O3 -ffast-math -msse2 -Ftree-vectorize -std = gnu99.

Für die Prüfung Ich lasse n eine Million sein, die Arrays so a[b[i]] Initialisierung und c[i] beide gleich 1.0/(i+1), mit einem paar verschiedenen Mustern von Indizes. Ein zuordnet a mit einer Million Elementen und setzt b in eine zufällige Permutation, eine andere zuordnet a mit 10M-Elemente und nutzt jeden 10., und die letzten zuordnet a mit 10M Elementen und setzt bis b[i+1] durch eine Zufallszahl von 1 bis 9 zu b[i] hinzufügen. Ich bin Timing, wie lange ein Anruf mit gettimeofday nimmt, die Caches Clearing clflush über den Arrays Aufruf und beträgt 1000 Versuche mit jeder Funktion. Ich geplottet Laufzeitverteilungen geglättet einige Codes aus den Eingeweiden von Kriterium (insbesondere die Kerndichteschätzer im statistics Paket).

Bandbreite

Nun, für die tatsächliche wichtige Anmerkung über die Bandbreite. 5333MB / s mit 2,4 GHz Takt ist etwas mehr als zwei Bytes pro Zyklus. Meine Daten ist lang genug, dass nichts zwischenspeicherbar sein sollte, und die Multiplikation der Laufzeit meiner Schleife durch (16 + 2 * 16 + 4 * 64) pro Iteration geladen Bytes, wenn alles vermisst mich gibt fast genau die ~ 5333MB / s Bandbreite mein System hat . Es sollte ziemlich einfach sein, dass die Bandbreite ohne SSE zu sättigen. Selbst unter der Annahme a wurden vollständig im Cache gespeichert, nur b und c für eine Iteration bewegt 12 Bytes von Daten zu lesen, und die naiven eine neue Iteration jemals dritten Zyklus mit Pipelining beginnen kann.

alles Unter der Annahme, weniger als die vollständige Caching auf a macht Arithmetik und Instruktion Zählung noch weniger einen Engpass. Ich wäre nicht überrascht, wenn die meisten der Speedup in meinem Code kommen von weniger Lasten b Ausgabe und c so mehr Platz frei zu verfolgen und zu spekulieren Vergangenheit Cache-Misses auf a.

Wider Hardware könnte mehr Unterschied machen. Ein Nehalem System drei Kanäle von DDR3-1333 laufen müßte 10667 * 3 / 2,66 = 12,6 Bytes pro Zyklus zu sättigen, die Speicherbandbreite bewegen. Das wäre für einen einzelnen Thread unmöglich sein, wenn a passt in Cache - aber bei 64 eine Linie Cache-Misses auf dem Vektor-Bytes summieren sich schnell - nur eine der vier Lasten in meiner Schleife in Caches fehlt die durchschnittliche erforderliche Bandbreite auf 16 Bytes bringt / Zyklus.

Answer 4

kurze Antwort nein. Lange Antwort ja, aber nicht effizient. Sie werden die Strafe dafür nicht ausgerichteten Belastungen entstehen, die jede Art von Vorteil zunichte machen wird. Es sei denn, Sie, dass b garantieren kann [i] aufeinanderfolgende Indizes ausgerichtet sind, werden Sie wahrscheinlich haben schlechtere Leistung nach der Vektorisierung

Wenn Sie vorher wissen, was Indizes sind, die beste, die zu entrollen und gibt expliziten Indizes. Ich habe etwas ähnliches Template-Spezialisierung und Code-Generierung verwendet wird. Wenn Sie interessiert sind, kann ich teilen

Ihr Kommentar zu beantworten, haben Sie grundsätzlich zu konzentrieren sich auf ein Array. Einfachste Sache sofort versuchen Sie Schleife um den Faktor zwei, Last niedrig und hoch ein separat zu blockieren, und dann mit mm * _ pd wie gewöhnlich. Pseudo-Code:

__m128d a, result;
for(i = 0; i < n; i +=2) {
  ((double*)(&a))[0] = A[B[i]];
  ((double*)(&a))[1] = A[B[i+1]];
  // you may also load B using packed integer instruction
  result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i])));
}

Ich erinnere mich nicht genau Funktionsnamen, zu überprüfen möchte. Außerdem beschränkt die Verwendung mit den Zeigern Schlüsselwort, wenn Sie es wissen kann keine Aliasing-Probleme sein. Dies ermöglicht es Compiler viel aggressiver zu sein.

Answer 5

Dies wird nicht vektorisiert werden, wie es ist, wegen der doppelten Indirektheit des Array-Indizes. Da Sie mit Doppel arbeiten sind es wenig oder nichts von SSE gewonnen werden, vor allem, da die meisten modernen CPUs haben 2 FPU trotzdem.