É possível vetorizar myNum += a[b[i]] * c[i];em x86_64?

Question

Quais intrínsecos eu usaria para vetorizar o seguinte (se é possível vetorizar) no 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

Aqui está minha tentativa, totalmente otimizada e testada:

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

Isso produz um código de montagem muito bonito usando gcc -O2 -msse2 (4.4.1).

Como você pode dizer, ter um uniforme n fará com que esse loop seja mais rápido e alinhado c. Se você pode se alinhar c, mudança _mm_loadu_pd para _mm_load_pd para um tempo de execução ainda mais rápido.

Answer 2

Eu começaria desenrolando o loop. Algo como

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;

Espero que isso permita ao compilador intercalar as cargas com a aritmética; perfil e observe a montagem para ver se há uma melhoria. Idealmente, o compilador gerará instruções SSE, mas não sou se isso acontecer na prática.

A desenrolar pode permitir que você faça isso:

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

(Desculpas pelo pseudocódigo no início e no fim, acho que a parte importante era o loop). Não sei ao certo se isso será mais rápido; Depende das várias latências e de quão bem o compilador pode reorganizar tudo. Certifique -se de perfilar antes e depois para ver se houve uma melhoria real.

Espero que ajude.

Answer 3

Os processadores Intel podem emitir duas operações de ponto flutuante, mas uma carga por ciclo, portanto, o acesso à memória é a restrição mais rígida.Com isso em mente, pretendi primeiro usar cargas compactadas para reduzir o número de instruções de carga e usei aritmética compactada apenas porque era conveniente.Desde então, percebi que a saturação da largura de banda da memória pode ser o maior problema, e toda a confusão com as instruções SSE pode ter sido uma otimização prematura se o objetivo fosse acelerar o código em vez de aprender a vetorizar.

SSE

O menor número de cargas possíveis sem nenhuma suposição sobre os índices em b requer desenrolar o loop quatro vezes.Uma carga de 128 bits obtém quatro índices de b, duas cargas de 128 bits obtêm cada uma um par de duplos adjacentes de c, e reunindo a exigia cargas independentes de 64 bits.Isso representa um piso de 7 ciclos por quatro iterações para código serial.(o suficiente para saturar a largura de banda da minha memória se o acesso a a não armazena em cache bem).Deixei de fora algumas coisas irritantes, como lidar com uma série de iterações que não são múltiplas de 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

Obter os índices é a parte mais complicada. movdqa carrega 128 bits de dados inteiros de um endereço alinhado de 16 bytes (Nehalem tem penalidades de latência para misturar as instruções SSE "inteiras" e "flutuantes"). punpckhqdq move 64 bits altos para 64 bits baixos, mas no modo inteiro, ao contrário do nome mais simples movhlpd.Mudanças de 32 bits são feitas nos registradores de uso geral. movhpd carrega um duplo na parte superior de um registro xmm sem perturbar a parte inferior - isso é usado para carregar os elementos de a diretamente em registradores compactados.

Este código é distintamente mais rápido que o código acima, que por sua vez é mais rápido que o código simples, e em todos os padrões de acesso, exceto no caso simples B[i] = i onde o loop ingênuo é realmente mais rápido.Eu também tentei algumas coisas como uma função em torno SUM(A(B(:)),C(:)) em Fortran que acabou basicamente equivalente ao loop simples.

Testei em um Q6600 (65 nm Core 2 a 2,4 Ghz) com 4 GB de memória DDR2-667, em 4 módulos.Testar a largura de banda da memória fornece cerca de 5.333 MB/s, então parece que estou vendo apenas um único canal.Estou compilando com o gcc 4.3.2-1.1 do Debian, -O3 -Ffast-math -msse2 -Ftree-vectorize -std=gnu99.

Para testar estou deixando n ser um milhão, inicializando os arrays para a[b[i]] e c[i] ambos iguais 1.0/(i+1), com alguns padrões diferentes de índices.Um aloca a com um milhão de elementos e conjuntos b para uma permutação aleatória, outro aloca a com 10 milhões de elementos e usa a cada 10, e o último aloca a com 10 milhões de elementos e configurações b[i+1] adicionando um número aleatório de 1 a 9 para b[i].Estou cronometrando quanto tempo leva uma ligação com gettimeofday, limpando os caches chamando clflush sobre as matrizes e medindo 1000 tentativas de cada função.Eu plotei distribuições de tempo de execução suavizadas usando algum código das entranhas de critério (em particular, o estimador de densidade do kernel no statistics pacote).

Largura de banda

Agora, para uma observação importante sobre largura de banda.5333 MB/s com clock de 2,4 GHz equivale a pouco mais de dois bytes por ciclo.Meus dados são longos o suficiente para que nada possa ser armazenado em cache, e multiplicar o tempo de execução do meu loop por (16 + 2 * 16 + 4 * 64) bytes carregados por iteração se tudo falhar me dá quase exatamente a largura de banda de ~ 5333 MB/s que meu sistema tem .Deve ser muito fácil saturar essa largura de banda sem SSE.Mesmo assumindo a foram completamente armazenados em cache, apenas lendo b e c para uma iteração move 12 bytes de dados, e o ingênuo pode iniciar uma nova iteração a cada terceiro ciclo com pipeline.

Assumindo algo menos do que o cache completo em a faz com que a aritmética e as instruções sejam ainda menos gargalos.Eu não ficaria surpreso se a maior parte da aceleração do meu código viesse da emissão de menos cargas para b e c portanto, há mais espaço livre para rastrear e especular falhas de cache anteriores a.

Hardware mais amplo pode fazer mais diferença.Um sistema Nehalem executando três canais de DDR3-1333 precisaria mover 3*10667/2,66 = 12,6 bytes por ciclo para saturar a largura de banda da memória.Isso seria impossível para um único thread se a cabe no cache - mas em 64 bytes um cache de linha falha no vetor aumenta rapidamente - apenas uma das quatro cargas em meu loop ausentes nos caches aumenta a largura de banda média necessária para 16 bytes/ciclo.

Answer 4

Resposta curta não. Resposta longa sim, mas não com eficiência. Você sofrerá a penalidade por fazer cargas não alinhadas, o que negará qualquer tipo de benefício. A menos que você possa garantir que os índices sucessivos b [i] estejam alinhados, você provavelmente terá um desempenho pior após a vetorização

Se você souber de antemão quais são os índices, o seu melhor é revelar e especificar índices explícitos. Fiz algo semelhante usando especialização de modelos e geração de código. Se você estiver interessado, eu posso compartilhar

Para responder ao seu comentário, você basicamente precisa se concentrar em uma matriz. A coisa mais fácil de tentar imediatamente é bloquear o loop de um fator de dois, carregar baixo e alto separadamente e depois usar milímetros*_pd como normalmente. Pseudo-código:

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

Não me lembro exatamente de nomes de funções, posso querer verificar novamente. Além disso, use a palavra -chave restrita com os ponteiros, se você souber que não pode haver problemas de alias. Isso permitirá que o compilador seja muito mais agressivo.

Answer 5

Isso não vai vetorizar como é, por causa da dupla indireção dos índices de matriz. Como você está trabalhando com duplas, há pouco ou nada a ser ganho da SSE, principalmente porque a maioria das CPUs modernas tem 2 FPUs de qualquer maneira.