Acelere a multiplicação da matriz 4x4 flutuante com SIMD

Muitos anos se passaram desde que me familiarizei com as instruções MMX, SSE e AVX posterior nos processadores Intel. Ao mesmo tempo, eles pareciam algum tipo de mágica contra o pano de fundo da montadora x86, que havia sido algo mundano. Eles me enganaram tanto que, alguns anos atrás, tive a idéia de escrever meu próprio renderizador de software para um jogo famoso. O que prometeu essas instruções me prometeu isso. Em algum momento, eu até pensei em escrever. Mas escrever texto acabou sendo muito mais complicado que código.

Naquela época, eu queria evitar problemas de suporte em diferentes processadores. Queria poder verificar meu representante na quantidade máxima disponível. Eu ainda tenho amigos com os antigos processadores AMD, e seu teto era SSE3. Portanto, naquele momento, decidi me limitar a um máximo de SSE3. Portanto, havia uma biblioteca matemática vetorial, um pouco menos do que totalmente implementada no SSE, com uma rara inclusão antes do SSE3. No entanto, em algum momento, pensei em qual desempenho máximo eu poderia extrair do processador para várias operações matemáticas vetoriais críticas. Uma dessas operações é multiplicar 4 matrizes flutuantes por 4.


Na verdade, eu decidi fazer esse negócio mais pelo bem do entretenimento. Eu já escrevi e tenho usado a multiplicação de matrizes para a renderização do meu software no SSE e isso parece ser suficiente para mim. Mas então decidi ver quantas medidas posso extrair em princípio da multiplicação de 2 matrizes float4x4. No meu SSE atual, são 16 ciclos de relógio. É verdade que a recente transição para a IACA 3 começou a mostrar 19, desde que comecei a escrever 1 * para algumas instruções em vez de 0 *. Aparentemente, antes, era apenas uma falha no analisador.

Brevemente sobre os utilitários usados

Para análise de código, usei o famoso utilitário Intel Architecture Code Analyzer . Para análise, uso a arquitetura Haswell (HSW), no mínimo, com suporte ao AVX2. Escrever código também é muito conveniente de usar: Guia Intel Intrinsics e manual de otimização da Intel .

Para montagem, uso a Comunidade MSVS 2017 no console. Eu escrevo o código na versão com intrínsecas. Você escreve uma vez e geralmente funciona imediatamente em plataformas diferentes. Além disso, o compilador x64 VC ++ não suporta assembler embutido, mas também quero que ele funcione no x64.

Como este artigo já está um pouco além do nível do iniciante na programação do SIMD, não descreverei registros, instruções, desenharei (ou faça a barba) belas imagens e tentarei aprender a programação usando as instruções do SIMD. O site da Intel está cheio de documentação excelente, clara e detalhada.

Eu queria tornar tudo mais fácil ... Mas acabou como sempre

É aqui que o momento começa, o que complica muito a implementação e o artigo. Portanto, vou me debruçar um pouco sobre isso. Não é interessante escrever multiplicação de matrizes com um layout de linha padrão de elementos. Quem precisava, e por isso estudaram nas universidades ou por conta própria. Nosso objetivo é a produtividade. Em primeiro lugar, mudei para o layout da coluna há muito tempo. Meu renderizador de software é baseado na API OpenGL e, portanto, para evitar transposições desnecessárias, comecei a armazenar elementos em colunas. Isso também é importante porque a multiplicação da matriz não é tão crítica. Matrizes 2-5-10 bem multiplicadas. E é isso. E então multiplicamos a matriz finalizada por milhares ou milhões de vértices. E esta operação é muito mais crítica. Você pode, é claro, transpor toda vez. Mas por que, se isso puder ser evitado.

Mas voltando às matrizes exclusivamente. Determinamos o armazenamento em colunas. No entanto, pode ser ainda mais complicado. É mais conveniente para mim armazenar os elementos seniores dos vetores e linhas da matriz nos registros SIMD para que x esteja no ponto mais alto (índice 3) ew no menor (índice 0). Aqui, aparentemente, teremos que recuar novamente sobre o porquê.

O problema é que, em um renderizador de software em um vetor, você precisa manipular o componente w com mais frequência ( 1 / z é armazenado lá), e é muito conveniente fazer isso através da versão _ss da operação (operações exclusivamente com o componente no flutuador inferior do registro xmm ), sem tocar em x, y, z . Portanto, no registro SSE, o vetor é armazenado em uma ordem compreensível x, y, z, w e na memória no inverso w, z, y, x .

Além disso, todas as opções de multiplicação também são implementadas por funções individuais. Isso é feito porque eu os uso para substituir a opção desejada, dependendo do tipo de instruções suportadas. Bem descrito aqui.

Implementamos a funcionalidade básica

Multiplicação com loops, linha ordenada

Tudo é simples e claro aqui. Para cada elemento, fazemos 4 multiplicações e 3 adições. No total, são 64 multiplicações e 48 adições. E isso sem levar em consideração a leitura dos elementos do registro.

Tudo é triste, em suma. Para esta opção, para o ciclo interno, a IACA emitiu: 3,65 ciclos de relógio para montagem x86 e 2,97 relógios para montagem x64 . Não pergunte por que números fracionários. Eu não sei A IACA 2.1 não sofreu com isso. De qualquer forma, esses números devem ser multiplicados por cerca de 4 * 4 * 4 = 64. Mesmo que você faça x64, o resultado será de 192 medidas. É claro que esta é uma estimativa aproximada. Não vejo o ponto de avaliar o desempenho com mais precisão para esta opção.

Implementação em loop, coluna ordenada

Multiplicação do ciclo, armazenamento orientado a SIMD

Essa implementação simplifica um pouco a compreensão do que está acontecendo por dentro, mas claramente não é suficiente.

Classes auxiliares

Para a conveniência de entender e escrever o código de referência e depuração, é conveniente implementar algumas classes auxiliares. Nada mais, tudo é apenas para entender. Observo que a implementação de classes de vetores e matrizes completas é uma questão difícil separada e não está incluída no tópico deste artigo.

Função de referência para testes

Como a ordem aceita dos elementos na matriz complica muito o entendimento, também não nos incomodaremos com a função compreensível de referência, que mostrará em futuras implementações que tudo funciona corretamente. Vamos comparar os resultados subsequentes com ele.


O algoritmo clássico é claramente pintado aqui, é difícil cometer um erro (mas você pode :-)). Nela, a IACA emitia: medidas x86 - 69,95, medidas x64 - 64 . Aqui estão cerca de 64 ciclos e veremos a aceleração dessa operação no futuro.

Implementação SSE

Algoritmo SSE clássico

Por que clássico? Porque há muito tempo está implementando o FVec como parte do MSVS. Para começar, escreveremos como apresentamos os elementos da matriz nos registros SSE. Já parece mais simples aqui. Apenas uma matriz transposta.

 //     00, 10, 20, 30 // m[0] -  SIMD /   01, 11, 21, 31 // m[1] 02, 12, 22, 32 // m[2] 03, 13, 23, 33 // m[3] 

Pegamos o código de desenrolamento da variante acima. De alguma forma, ele é hostil à SSE. O primeiro grupo de linhas consiste nos resultados da coluna da matriz resultante: r._00, r._01, r._02, r._03 . Temos essa coluna, mas precisamos de uma linha. Sim, m , n parecem inconvenientes para os cálculos. Portanto, reorganizamos as linhas do algoritmo para que o resultado r seja de linha.

 //  ,     r[0] r00 = m00*n00 + m01*n10 + m02*n20 + m03*n30; r10 = m10*n00 + m11*n10 + m12*n20 + m13*n30; r20 = m20*n00 + m21*n10 + m22*n20 + m23*n30; r30 = m30*n00 + m31*n10 + m32*n20 + m33*n30; //  ,     r[1] r01 = m00*n01 + m01*n11 + m02*n21 + m03*n31; r11 = m10*n01 + m11*n11 + m12*n21 + m13*n31; r21 = m20*n01 + m21*n11 + m22*n21 + m23*n31; r31 = m30*n01 + m31*n11 + m32*n21 + m33*n31; //  ,     r[2] r02 = m00*n02 + m01*n12 + m02*n22 + m03*n32; r12 = m10*n02 + m11*n12 + m12*n22 + m13*n32; r22 = m20*n02 + m21*n12 + m22*n22 + m23*n32; r32 = m30*n02 + m31*n12 + m32*n22 + m33*n32; //  ,     r[3] r03 = m00*n03 + m01*n13 + m02*n23 + m03*n33; r13 = m10*n03 + m11*n13 + m12*n23 + m13*n33; r23 = m20*n03 + m21*n13 + m22*n23 + m23*n33; r33 = m30*n03 + m31*n13 + m32*n23 + m33*n33; 

Mas isso já é muito melhor. O que, de fato, vemos? De acordo com as colunas do algoritmo em cada grupo, temos as linhas da matriz m envolvidas:

 m [0] = {00,10,20,30}, m [1] = {01,11,21,31}, m [2] = {02,12,22,32}, m [3] = {03,13,23,33},

que são multiplicados pelo mesmo elemento da matriz n . Por exemplo, para o primeiro grupo é: n._00, n._10, n._20, n._30 . E os elementos da matriz n para cada grupo de linhas do algoritmo estão novamente em uma linha da matriz.

Então, tudo é simples: simplesmente pegamos as linhas da matriz m por índice, mas, quanto aos elementos n , pegamos sua linha e a embaralhamos pelos 4 elementos do registro através da instrução shuffle , para multiplicar pela linha da matriz m no registro. Por exemplo, para o elemento n._00 (lembre-se de que sua mudança no registro possui o índice 3), será:

  _mm_shuffle_ps (n [0], n [0], _MM_SHUFFLE (3,3,3,3)) 

De uma forma simplificada, o algoritmo se parece com isso:

 //   n[0]={00,10,20,30} r[0] = m[0] * n00 + m[1] * n10 + m[2] * n20 + m[3] * n30; //   n[1]={01,11,21,31} r[1] = m[0] * n01 + m[1] * n11 + m[2] * n21 + m[3] * n31; //   n[2]={02,12,22,32} r[2] = m[0] * n02 + m[1] * n12 + m[2] * n22 + m[3] * n32; //   n[3]={03,13,23,33} r[3] = m[0] * n03 + m[1] * n13 + m[2] * n23 + m[3] * n33; 

Agora, alteramos os elementos n no algoritmo para o shuffle correspondente, multiplicação por _mm_mul_ps , a soma por _mm_add_ps , e pronto. Isso funciona. O código parece, no entanto, muito pior do que o próprio algoritmo parecia. Nesse código, a IACA emitia: x86 - 18,89, x64 - 16 ciclos . Isso é 4 vezes mais rápido que o anterior. No SSE registre o 4º carro alegórico Relação quase linear.

Decore a implementação do SSE

Ainda assim, no código, parece horrível. Vamos tentar melhorar isso escrevendo um pouco de açúcar sintático.


O compilador pode alinhar perfeitamente essas funções (embora algumas vezes sem o __forceinline de qualquer forma).


E, portanto, já é muito melhor e mais legível. Para isso, a IACA produziu aproximadamente o resultado esperado: x86 - 19 (e por que não fracionário?), X64 - 16 . De fato, o desempenho não mudou, mas o código é muito mais bonito e compreensível.

Pouca contribuição para otimização futura

Vamos introduzir mais uma melhoria no nível de uma função que apareceu recentemente na versão de ferro. Operação de adição múltipla (fma) . fma (a, b, c) = a * b + c .


Por que isso é necessário? Primeiro de tudo, para otimização futura. Por exemplo, você pode simplesmente substituir mad por fma no código finalizado através das mesmas macros que desejar. Mas vamos estabelecer as bases para a otimização agora:


IACA: x86 - 18,89, x64 - 16 . Novamente fracionário. Ainda assim, a IACA às vezes produz resultados estranhos. O código não mudou muito. Provavelmente até um pouco pior. Mas a otimização às vezes exige tais sacrifícios.

Passamos à economia através de _mm_stream

Vários guias de otimização recomendam mais uma vez não puxar o cache para operações de salvamento em massa. Isso geralmente é justificado quando você está processando vértices de milhares ou mais. Mas para matrizes, isso talvez não seja tão importante. No entanto, vou adicioná-lo de qualquer maneira.


Nada mudou com o tempo aqui, a partir da palavra. Mas, de acordo com as recomendações, agora não tocamos no cache mais uma vez.

Implementação AVX

Opção AVX básica

Em seguida, passamos ao próximo estágio de otimização. O quarto float está incluído no registro SSE e no AVX já é 8. Ou seja, existe uma chance teórica de reduzir o número de operações executadas e aumentar a produtividade, se não pela metade, pelo menos 1,5 vezes. Mas algo me diz que nem tudo será tão simples com a transição para o AVX. Podemos obter os dados necessários de registros duplos?

Vamos tentar descobrir. Novamente, escrevemos nosso algoritmo de multiplicação usado acima. Você não pode fazer isso, mas é mais conveniente lidar com o código quando tudo estiver próximo e você não precisará rolar meia página.

 //    : 00, 10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33 //   SSE: r0 = m0*n00 + m1*n10 + m2*n20 + m3*n30 r1 = m0*n01 + m1*n11 + m2*n21 + m3*n31 r2 = m0*n02 + m1*n12 + m2*n22 + m3*n32 r3 = m0*n03 + m1*n13 + m2*n23 + m3*n33 

Na saída, esperamos obter o resultado em ymm = {r0: r1} e ymm = {r2: r3} . Se na versão SSE nosso algoritmo foi generalizado para colunas, agora precisamos generalizá-lo para linhas. Portanto, agir como no caso da opção SSE não funcionará.

Se considerarmos a matriz m nos registradores ymm , obtemos ymm = {m0: m1} e ymm = {m2: m3}, respectivamente. Anteriormente, tínhamos apenas colunas matriciais no registro e agora colunas e linhas.

Se você tentar agir como antes, multiplique ymm = {m0: m1} pelo registro ymm = {n00, n00, n00, n00}: {n10, n10, n10, n10} . Como n00 e n01 estão na mesma linha da matriz n , a julgar pelo conjunto disponível de instruções AVX, espalhá-las por ymm será caro. O shuffle e o permute funcionam separadamente para cada um dos quatro quatros de um flutuador (alto e baixo xmm ) nos registros de ymm .

Se tirarmos ymm da matriz n , obteremos os dois elementos n00 e n10 no máximo de 2 xmm dentro do registro ymm . {n00, n10, n20, n30}: {n01, n11, n21, n31} . Geralmente, o índice para as instruções existentes é de 0 a 3. E ele aborda float apenas dentro de um registro xmm dentre dois registros ymm internos. Não é possível transferir o n10 do xmm mais antigo para o mais novo mais barato . E então esse foco deve ser repetido várias vezes. Não podemos suportar essa perda de medidas. É necessário inventar outra coisa.

Costumávamos generalizar colunas, mas agora linhas. Portanto, tentaremos ir de uma maneira um pouco diferente. Precisamos obter o resultado em {r0: r1} . Isso significa que o algoritmo deve ser aprimorado não em linhas separadas do algoritmo, mas em duas de uma vez. E aqui, o que foi menos no trabalho de embaralhar e permutar , será uma vantagem para nós. Observamos o que teremos nos registros ymm quando considerarmos a matriz n .

 n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31} n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33} 

Sim, percebemos que em diferentes partes xmm do registro ymm temos os elementos 00 e 01 . Eles podem ser multiplicados caso por registro através do comando permute em {_00, _00, _00, _00}: {_ 01, _01, _01, _01} , indicando apenas um índice 3 para ambas as partes xmm . É exatamente disso que precisamos. De fato, os coeficientes também são usados ​​em diferentes linhas. Somente agora no registro ymm correspondente para multiplicação será necessário manter {m0: m0} , ou seja, a primeira linha duplicada da matriz m .

Então, pintamos o algoritmo com mais detalhes. Lemos as linhas duplas da matriz m em registros ymm :

 mm[0] = {m0:m0} mm[1] = {m1:m1} mm[2] = {m2:m2} mm[3] = {m3:m3} 

E então calcularemos a multiplicação como:

 r0r1 = mm[0] * {n00,n00,n00,n00:n01,n01,n01,n01} + // permute<3,3,3,3>(n0n1) mm[1] * {n10,n10,n10,n10:n11,n11,n11,n11} + // permute<2,2,2,2>(n0n1) mm[2] * {n20,n20,n20,n20:n21,n21,n21,n21} + // permute<1,1,1,1>(n0n1) mm[3] * {n30,n30,n30,n30:n31,n31,n31,n31} // permute<0,0,0,0>(n0n1) r2r3 = mm[0] * {n02,n02,n02,n02:n03,n03,n03,n03} + // permute<3,3,3,3>(n2n3) mm[1] * {n12,n12,n12,n12:n13,n13,n13,n13} + // permute<2,2,2,2>(n2n3) mm[2] * {n22,n22,n22,n22:n23,n23,n23,n23} + // permute<1,1,1,1>(n2n3) mm[3] * {n32,n32,n32,n32:n33,n33,n33,n33} // permute<0,0,0,0>(n2n3) 

Reescrevemos mais claramente:

 r0r1 = mm[0]*n0n1<3,3,3,3>+mm[1]*n0n1<2,2,2,2>+mm[2]*n0n1<1,1,1,1>+mm[3]*n0n1<0,0,0,0> r2r3 = mm[0]*n2n3<3,3,3,3>+mm[1]*n2n3<2,2,2,2>+mm[2]*n2n3<1,1,1,1>+mm[3]*n2n3<0,0,0,0> 

Ou de forma simplificada:

 r0r1 = mm[0]*n0n1<3> + mm[1]*n0n1<2> + mm[2]*n0n1<1> + mm[3]*n0n1<0> r2r3 = mm[0]*n2n3<3> + mm[1]*n2n3<2> + mm[2]*n2n3<1> + mm[3]*n2n3<0> 

Tudo parece estar claro.


Aqui estão os números interessantes da IACA: x86 - 12,53, x64 - 12 . Embora, é claro, eu quisesse melhor. Perdeu alguma coisa.

Otimização AVX mais açúcar sintático

Parece que no código acima, o AVX não estava acostumado a todo o seu potencial. Descobrimos que, em vez de definir duas linhas idênticas no registro ymm , podemos usar broadcast , que pode preencher o registro ymm com dois valores xmm idênticos. Além disso, ao longo do caminho, adicione um pouco de "açúcar sintático" para as funções do AVX.


E aqui os resultados já são mais interessantes. A IACA produz números: x86 - 10, x64 - 8,58 , que parecem muito melhores, mas ainda não 2 vezes.

Opção AVX + FMA (final)

Vamos fazer mais uma tentativa. Agora, seria lógico recuperar o conjunto de instruções FMA novamente, pois foi adicionado aos processadores após o AVX. Basta alterar mul + add individual para uma operação. Embora ainda utilizemos a instrução de multiplicação para fornecer ao compilador mais oportunidades de otimização e ao processador para execução paralela de multiplicações. Normalmente, eu olho para o código gerado no assembler para garantir que opção é melhor.

Nesse caso, precisamos calcular a * b + c * d + e * f + g * h . Você pode fazer isso na testa: fma (a, b, fma (c, d, fma (e, f, g * h))) . Mas, como vemos, é impossível executar uma operação aqui sem concluir a anterior. E isso significa que não poderemos usar a capacidade de fazer multiplicações emparelhadas, como o pipeline do SIMD nos permite fazer. Se transformarmos ligeiramente os cálculos fma (a, b, c * d) + fma (e, f, g * h) , veremos que podemos paralelizar os cálculos. Primeiro faça duas multiplicações independentes e depois duas operações fma independentes.


IACA: x86 - 9.21, x64 - 8 . Agora está muito bom. Alguém provavelmente dirá o que pode ser feito ainda melhor, mas não sei como.

Benchmarks

Observo imediatamente que esses números não devem ser tomados como a verdade suprema. Mesmo com um teste fixo, eles nadam dentro de certos limites.E ainda mais, eles se comportam de maneira diferente em plataformas diferentes. Com qualquer otimização, faça medições especificamente para o seu caso.

Sumário

• Função: nome da função. O final em s - funciona com uma transmissão, mov diferentemente normal (sem transmissão). Adicionado para maior clareza, pois isso é importante o suficiente.
• Ciclos da IACA: número de ticks por função calculado pela IACA
• Ciclos medidos: número medido de medidas (menos é mais)
• Aceleração da IACA: número de medidas em uma linha zero / número de medidas em uma linha
• Aceleração medida: número de medidas na linha zero / número de medidas na linha (quanto mais, melhor)

Para loop_m, os ticks do artigo foram multiplicados por 64. Ou seja, esse é um valor muito aproximado. De fato, aconteceu dessa maneira.

i3-3770:

i7-8700K:

Código de teste na fonte. Se houver sugestões razoáveis ​​sobre como melhorá-las, escreva nos comentários.

BÔNUS do reino da ficção

Na verdade, é do reino da ficção, porque se eu vi processadores com suporte para AVX512, talvez nas fotos. No entanto, tentei implementar o algoritmo. Aqui não vou explicar nada, uma analogia completa com o AVX + FMA. O algoritmo é o mesmo, apenas menos operações.


Os números são fantásticos: x86 - 4,79, x64 - 5,42 (IACA com arquitetura SKX). Isso apesar do algoritmo ter 64 multiplicações e 48 adições.

Código PS do artigo

Esta é a minha primeira experiência escrevendo um artigo. Obrigado a todos por seus comentários. Eles ajudam a melhorar o código e o artigo.