E 'possibile vettorizzare myNum + = a [b [i]] * c [i]; su x86_64?

Question

Cosa intrinseche dovrei usare per vettorizzare il seguente (se è ancora possibile vettorizzare) sul 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

Ecco il mio andare a esso, completamente ottimizzato e testato:

#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))
));

Questo produce molto bella codice assembly utilizzando gcc -O2 -msse2 (4.4.1).

Come si può dire, con un ancora n farà questo circuito più veloce così come un c allineata. Se è possibile allineare c, cambiare _mm_loadu_pd per _mm_load_pd per un tempi di esecuzione più veloci.

Answer 2

Vorrei iniziare srotolando il ciclo. Qualcosa di simile

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;

Si spera che permette al compilatore di interleave i carichi con l'aritmetica; profilo e guardare l'assemblea a vedere se c'è un miglioramento. Idealmente il compilatore genera le istruzioni SSE, ma non sono se questo accade nella pratica.

apertolo ulteriormente potrebbe permetterà di fare questo:

__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...

(scuse per lo pseudocodice all'inizio e alla fine, ho dato la parte importante è stata la loop). Non so con certezza se quello sarà più veloce; dipende dalle varie latenze e quanto bene il compilatore può riorganizzare tutto. Assicurati di profilo prima e dopo per vedere se v'è stato un miglioramento effettivo.

La speranza che aiuta.

Answer 3

processori Intel possono emettere due operazioni in virgola mobile ma un carico per ciclo, così accessi alla memoria sono il vincolo stretto. Con questo in mente, ho puntato primo ad utilizzare carichi imballati per ridurre il numero delle istruzioni di carico, e usato confezionato aritmetica solo perché era conveniente. Da allora ho capito che saturare la larghezza di banda della memoria può essere il più grande problema, e tutto il fare in giro con istruzioni SSE avrebbe potuto essere l'ottimizzazione prematura se il punto è stato quello di rendere il codice andare veloce, piuttosto che imparare a vettorizzare.

SSE

il minor numero di possibili carichi senza assunzione sugli indici in b richiede srotolando il ciclo quattro volte. Un carico 128 bit ottiene quattro indici da b, due carichi 128 bit ciascuno ottenere un accoppiamento dei doppi adiacenti da c, e la raccolta a richiesto carichi indipendenti 64 bit. Questo è un pavimento di 7 cicli per quattro iterazioni per codice seriale. (Abbastanza per saturare la mia larghezza di banda di memoria se l'accesso a a non memorizza nella cache bene). Ho lasciato fuori alcune cose fastidiose come la gestione di un numero di iterazioni che non è un multiplo di 4.

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

Come gli indici fuori è la parte più complicata. carichi movdqa 128 bit di dati interi da un indirizzo allineato 16 byte (Nehalem trovi sanzioni latenza per la miscelazione del "integer" e "float" istruzioni SSE). punpckhqdq muove alti 64 bit a 64 bit bassi, ma in modalità ad interi differenza del movhlpd più denominata semplicemente. 32 turni di bit sono fatte nei registri di uso generale. carichi movhpd una doppia nella parte superiore di un registro xmm senza disturbare la parte inferiore -. Questo è utilizzato per caricare elementi di a direttamente nei registri confezionati

Questo codice nettamente più veloce rispetto al codice di cui sopra, che è a sua volta più veloce del codice semplice, e su ogni modello di accesso, ma il semplice caso in cui il ciclo B[i] = i ingenuo è in realtà più veloce. Ho anche provato una cosa pochi, come una funzione intorno SUM(A(B(:)),C(:)) in Fortran che ha finito sostanzialmente equivalente al semplice ciclo.

ho testato su un Q6600 (65 nm Core 2 a 2,4 GHz) con 4GB di memoria DDR2-667, in 4 moduli. Test banda di memoria dà circa 5333 MB / s, così sembra come sto solo vedendo un singolo canale. Sto compilando con gcc di Debian 4.3.2-1.1, -O3 -ffast-math -msse2 -Ftree-vectorize -std = gnu99.

Per la prova Sto lasciando n essere un milione, l'inizializzazione degli array in modo a[b[i]] e c[i] sia uguale 1.0/(i+1), con una serie di vari modelli di indici. Una alloca a con un milione di elementi e imposta b ad una permutazione casuale, un altro alloca a con elementi 10M e utilizza ogni 10, e gli ultimi alloca a con elementi 10M e imposta b[i+1] aggiungendo un numero casuale da 1 a 9 per b[i]. Sono tempi quanto tempo una chiamata porta con gettimeofday, aprendo le cache chiamando clflush nel corso degli array, e, di 1000 prove di ogni funzione. Ho tracciato lisciai distribuzioni di runtime utilizzando codice dalle viscere di criterio di (in particolare, lo stimatore densità di kernel nella confezione statistics).

banda

Ora, per l'importante nota reale circa la larghezza di banda. 5333MB / s con orologio 2.4Ghz è poco più di due byte per ciclo. I miei dati è abbastanza lungo che nulla deve essere memorizzabile nella cache, e moltiplicando il tempo di esecuzione del mio loop (16 + 2 * 16 + 4 * 64) byte caricati per ogni iterazione se tutto manca mi dà quasi esattamente la larghezza di banda ~ 5333MB / s il mio sistema ha . Dovrebbe essere abbastanza facile per saturare la larghezza di banda senza SSE. Anche supponendo a erano completamente nella cache, semplicemente leggendo b e c per un'iterazione muove 12 byte di dati, e il naif può iniziare una nuova iterazione terzo ciclo mai con pipelining.

Supponendo niente di meno che la cache completa su a rende l'aritmetica e l'istruzione conta ancora meno di un collo di bottiglia. Non sarei sorpreso se la maggior parte del aumento di velocità nel mio codice viene dal rilascio di un minor numero di carichi da b e c così più spazio è libero di monitorare e speculare cache miss passati su a.

hardware più ampia potrebbe avere più differenza. Un sistema Nehalem esecuzione tre canali di DDR3-1333 dovrebbe muoversi 3 * 10667 / 2,66 = 12,6 byte per ciclo di saturare banda di memoria. Sarebbe impossibile per un singolo thread, se a si inserisce nella cache - ma a 64 byte di una cache linea manca sul vettore aggiungere rapidamente - solo uno dei quattro carichi nel mio anello mancante nella cache porta in primo piano la larghezza di banda media richiesta a 16 byte / ciclo.

Answer 4

risposta breve no. risposta lunga sì, ma non in modo efficiente. Si incorrerà nella penalità per fare i carichi non allineati che negare qualsiasi tipo di beneficio. A meno che non si può garantire che b [i] indici successivi sono allineati, è molto probabile che avere prestazioni peggiori dopo la vettorializzazione

se si sa in anticipo ciò che gli indici sono, il tuo migliore che è srotolare e specificare gli indici espliciti. Ho fatto qualcosa di simile utilizzando il modello di specializzazione e la generazione del codice. se siete interessati, posso condividere

di rispondere il tuo commento, che, fondamentalmente, di concentrarsi su un array. cosa più facile da provare subito è quello di bloccare si esegue un ciclo di un fattore due, caricare bassa e alta una parte, e quindi utilizzare mm * _ pd come di solito. Pseudocodice:

__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])));
}

Non mi ricordo i nomi delle funzioni esattamente, può decidere di doppio controllo. Inoltre, utilizzare limitare parola chiave con i puntatori se si sa non ci possono essere problemi di aliasing. Questo permetterà compilatore di essere molto più aggressivo.

Answer 5

Questo non sta andando vectorize come è, a causa del doppio riferimento indiretto degli indici di matrice. Dal momento che si sta lavorando con i doppi c'è poco o nulla da guadagnare da SSE, in particolare per quanto la maggior parte delle moderne CPU hanno 2 FPU in ogni caso.