Vectorizing Modular Arithmetic

Question

Sto provando a scrivere codice di aggiunta vettoriale ragionevolmente veloce componente-saggio. Sto lavorando con (firmato, credo) interi a 64 bit.

La funzione è

void addRq (int64_t* a, const int64_t* b, const int32_t dim, const int64_t q) { 
    for(int i = 0; i < dim; i++) { 
     a[i] = (a[i]+b[i])%q; // LINE1 
    } 
} 

Sto compilando con icc -std=gnu99 -O3 (ICC in modo da poter usare SVML in seguito) su un Ivybridge (SSE4.2 e AVX, ma non AVX2).

La mia linea di base sta rimuovendo %q da LINE1. 100 chiamate di funzione (iterata) con dim=11221184 richiedono 1,6 secondi. ICC auto-vettorizza il codice per SSE; grande.

Voglio davvero fare aggiunte modulari però. Con %q, ICC non auto-vettorializza il codice e viene eseguito in 11,8 secondi (!). Anche ignorando l'auto-vettorizzazione per il precedente tentativo, questo sembra ancora eccessivo.

Poiché non ho AVX2, la vettorizzazione con SSE richiede SVML, che è forse il motivo per cui ICC non ha auto-vettorializzato. In ogni caso, ecco il mio tentativo di vettorizzare il ciclo interno:

__m128i qs = _mm_set1_epi64x(q); 
for(int i = 0; i < dim; i+=2) { 
    __m128i xs = _mm_load_si128((const __m128i*)(a+i)); 
    __m128i ys = _mm_load_si128((const __m128i*)(b+i)); 
    __m128i zs = _mm_add_epi64(xs,ys); 
    zs = _mm_rem_epi64(zs,qs); 
    _mm_store_si128((__m128i*)(a+i),zs); 
} 

Assemblea per il ciclo principale è:

..B3.4:       # Preds ..B3.2 ..B3.12 
    movdqa (%r12,%r15,8), %xmm0       #59.22 
    movdqa %xmm8, %xmm1         #60.14 
    paddq  (%r14,%r15,8), %xmm0       #59.22 
    call  __svml_i64rem2        #61.9 
    movdqa %xmm0, (%r12,%r15,8)       #61.36 
    addq  $2, %r15          #56.30 
    cmpq  %r13, %r15         #56.24 
    jl  ..B3.4  # Prob 82%      #56.24 

Così il codice è sempre Vettorializzare come previsto. So che potrei non ottenere un aumento di velocità 2x grazie a SVML, ma il codice funziona in 12,5 secondi, più lentamente che senza vettorizzazione! È davvero il meglio che si può fare qui?

Answer 1

Né SSE2 né AVX2 hanno istruzioni di divisione intera. Intel non è in grado di chiamare le funzioni intrinseche SVML poiché molte di esse sono funzioni complicate che si associano a diverse istruzioni e non solo a poche.

C'è un modo per fare una divisione più veloce (e modulo) con SSE2 o AVX2. Vedi questo documento Improved division by invariant integers. Fondamentalmente si precomputa un divisore e poi si fa la moltiplicazione. Precaricare il divisore richiede tempo ma per un certo valore di dim nel tuo codice dovrebbe vincere. Ho descritto questo metodo in modo più dettagliato qui SSE integer division? Ho anche implementato con successo questo metodo in un certo numero cercatore primo Finding lists of prime numbers with SIMD - SSE/AVX

Agner Fog implementa a 32 bit (ma non a 64-bit) divisione nel suo Vector Class utilizzando il metodo descritto in questo documento . Sarebbe un buon punto di partenza se vuoi un po 'di codice ma dovrai estenderlo a 64-bit.

Modifica: in base ai commenti di Mysticial e supponendo che gli input siano già ridotti, ho prodotto una versione per SSE. Se questo è compilato in MSVC, è necessario che sia in modalità 64 bit in quanto la modalità a 32 bit non supporta _mm_set1_epi64x. Questo può essere risolto per la modalità a 32 bit, ma non voglio farlo.

#ifdef _MSC_VER 
#include <intrin.h> 
#endif 
#include <nmmintrin.h>     // SSE4.2 
#include <stdint.h> 
#include <stdio.h> 

void addRq_SSE(int64_t* a, const int64_t* b, const int32_t dim, const int64_t q) { 
    __m128i q2 = _mm_set1_epi64x(q); 
    __m128i t2 = _mm_sub_epi64(q2,_mm_set1_epi64x(1)); 
    for(int i = 0; i < dim; i+=2) { 
     __m128i a2 = _mm_loadu_si128((__m128i*)&a[i]); 
     __m128i b2 = _mm_loadu_si128((__m128i*)&b[i]); 
     __m128i c2 = _mm_add_epi64(a2,b2); 
     __m128i cmp = _mm_cmpgt_epi64(c2, t2); 
     c2 = _mm_sub_epi64(c2, _mm_and_si128(q2,cmp)); 
     _mm_storeu_si128((__m128i*)&a[i], c2); 
    } 
} 

int main() { 
    const int64_t dim = 20; 
    int64_t a[dim]; 
    int64_t b[dim]; 
    int64_t q = 10; 

    for(int i=0; i<dim; i++) { 
     a[i] = i%q; b[i] = i%q; 
    } 
    addRq_SSE(a, b, dim, q); 
    for(int i=0; i<dim; i++) { 
     printf("%d\n", a[i]); 
    } 
}