Acelere la multiplicación de matriz flotante 4x4 con SIMD

Han pasado muchos años desde que me familiaricé con las instrucciones MMX, SSE y más tarde AVX en los procesadores Intel. En un momento, parecían una especie de magia en el contexto del ensamblador x86, que durante mucho tiempo había sido algo mundano. Me engancharon tanto que hace un par de años tuve la idea de escribir mi propio procesador de software para un juego famoso. Lo que prometió estas instrucciones me prometió esto. En algún momento, incluso pensé en escribirlo. Pero escribir texto resultó ser mucho más complicado que el código.

En ese momento, quería evitar problemas con el soporte en diferentes procesadores. Quería poder verificar mi procesador en la cantidad máxima disponible. Todavía tengo amigos con viejos procesadores AMD, y su techo era SSE3. Por lo tanto, en ese momento decidí limitarme a un máximo de SSE3. Así que había una biblioteca matemática vectorial, un poco menos que totalmente implementada en SSE, con una rara inclusión antes de SSE3. Sin embargo, en algún momento me pregunté qué rendimiento máximo podría obtener del procesador para una serie de operaciones matemáticas de vectores críticos. Una de estas operaciones es multiplicar flotador 4 matrices por 4.


En realidad, decidí hacer este negocio más por el entretenimiento. Ya he escrito y he estado usando la multiplicación de matrices para la representación de mi software en SSE y parece ser suficiente para mí. Pero luego decidí ver cuántas medidas puedo extraer en principio de multiplicar 2 matrices float4x4. En mi SSE actual, estos son 16 ciclos de reloj. Es cierto que la transición reciente a IACA 3 comenzó a mostrar 19, ya que comencé a escribir 1 * para algunas instrucciones en lugar de 0 *. Aparentemente antes era solo una falla en el analizador.

Brevemente sobre las utilidades utilizadas

Para el análisis de código utilicé la famosa utilidad Intel Architecture Code Analyzer . Para el análisis, uso la arquitectura Haswell (HSW), como mínimo con soporte para AVX2. Escribir código también es muy conveniente de usar: Intel Intrinsics Guide y Intel optimization manual .

Para el ensamblaje, uso MSVS 2017 Community desde la consola. Escribo el código en la versión con intrínsecos. Escribe una vez, y generalmente funciona de inmediato en diferentes plataformas. Además, el compilador x64 VC ++ no es compatible con el ensamblador en línea, pero también quiero que funcione con x64.

Dado que este artículo ya está un poco más allá del nivel de principiante en la programación SIMD, no describiré registros, instrucciones, dibujaré (o afeitaré) bellas imágenes e intentaré aprender a programar usando las instrucciones SIMD. El sitio web de Intel está lleno de documentación excelente, clara y detallada.

Quería hacer todo más fácil ... Pero resultó como siempre

Aquí es donde comienza el momento, lo que complica mucho tanto la implementación como el artículo. Por lo tanto, me detendré un poco. No es interesante para mí escribir multiplicación matricial con un diseño de fila estándar de elementos. Quién lo necesitaba, por lo que estudiaron en universidades o por su cuenta. Nuestro objetivo es la productividad. En primer lugar, cambié al diseño de columna hace mucho tiempo. Mi renderizador de software se basa en la API de OpenGL y, por lo tanto, para evitar transposiciones innecesarias, comencé a almacenar elementos en columnas. Esto también es importante porque la multiplicación de matrices no es tan crítica. Bien multiplicado 2-5-10 matrices. Y eso es todo. Y luego multiplicamos la matriz terminada por miles o millones de vértices. Y esta operación es mucho más crítica. Puedes, por supuesto, transponer cada vez. Pero por qué, si esto se puede evitar.

Pero volvamos a las matrices exclusivamente. Determinamos el almacenamiento en columnas. Sin embargo, puede ser aún más complicado. Es más conveniente para mí almacenar los elementos superiores de vectores y filas de matriz en registros SIMD para que x esté en el flotante más alto (índice 3) yw en el menor (índice 0). Aquí, aparentemente, tendremos que retirarnos nuevamente sobre por qué.

El problema es que en un renderizador de software en un vector tienes que manipular el componente w con más frecuencia ( 1 / z se almacena allí), y es muy conveniente hacerlo a través de la versión _ss de la operación (operaciones exclusivamente con el componente en el flotador inferior del registro xmm ), sin tocar x, y, z . Por lo tanto, en el registro SSE, el vector se almacena en un orden comprensible x, y, z, w , y en la memoria en el reverso w, z, y, x .

Además, todas las opciones de multiplicación también se implementan mediante funciones individuales. Esto se hace porque los uso para sustituir la opción deseada según el tipo de instrucciones admitidas. Bien descrito aquí.

Implementamos la funcionalidad básica.

Multiplicación con bucles, fila ordenada

Aquí todo es simple y claro. Para cada elemento hacemos 4 multiplicaciones y 3 sumas. En total, son 64 multiplicaciones y 48 adiciones. Y esto sin tener en cuenta la lectura de los elementos de registro.

Todo es triste, en resumen. Para esta opción, para el ciclo interno, IACA emitió: 3.65 ciclos de reloj para ensamblaje x86 y 2.97 relojes para ensamblaje x64 . No preguntes por qué números fraccionarios. No lo se IACA 2.1 no sufrió esto. En cualquier caso, estos números deben multiplicarse por aproximadamente 4 * 4 * 4 = 64. Incluso si toma x64, el resultado es de aproximadamente 192 medidas. Está claro que esta es una estimación aproximada. No veo el punto de evaluar el rendimiento con mayor precisión para esta opción.

Implementación de bucle, columna ordenada

Multiplicación de ciclos, almacenamiento orientado a SIMD

Esta implementación simplifica un poco la comprensión de lo que está sucediendo en el interior, pero claramente no es suficiente.

Clases de ayuda

Para la conveniencia de comprender y escribir referencias y códigos de depuración, es conveniente implementar un par de clases auxiliares. Nada más, todo es solo para comprender. Observo que la implementación de clases de vectores y matrices completas es una pregunta difícil por separado, y no se incluye en el tema de este artículo.

Función de referencia para pruebas

Dado que el orden aceptado de los elementos en la matriz complica mucho la comprensión, tampoco nos molestará la función comprensible de referencia, que mostrará en implementaciones futuras que todo funciona correctamente. Compararemos los resultados posteriores con él.


El algoritmo clásico está claramente pintado aquí, es difícil cometer un error (pero puedes :-)). En él, IACA emitió: x86 - 69.95 medidas, x64 - 64 medidas . Aquí hay unos 64 ciclos y veremos la aceleración de esta operación en el futuro.

Implementación de SSE

Algoritmo SSE clásico

¿Por qué clásico? Porque ha estado en la implementación de FVec como parte de MSVS. Para comenzar, escribiremos cómo presentamos los elementos de la matriz en los registros SSE. Ya parece más simple aquí. Solo una matriz transpuesta.

 //     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] 

Tomamos el código de desenrollado de la variante anterior. De alguna manera él es hostil para SSE. El primer grupo de filas consiste en los resultados para la columna de la matriz resultante: r._00, r._01, r._02, r._03 . Tenemos esta columna, pero necesitamos una fila. Sí, y m , n parecen inconvenientes para los cálculos. Por lo tanto, reorganizamos las líneas del algoritmo para que el resultado r sea fila-sabio.

 //  ,     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; 

Pero esto ya es mucho mejor. ¿Qué, de hecho, vemos? De acuerdo con las columnas del algoritmo en cada grupo, tenemos las filas de la matriz m involucradas:

 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 se multiplican por el mismo elemento de la matriz n . Por ejemplo, para el primer grupo es: n._00, n._10, n._20, n._30 . Y los elementos de la matriz n para cada grupo de filas del algoritmo nuevamente se encuentran en una fila de la matriz.

Entonces todo es simple: simplemente tomamos las filas de la matriz m por índice, pero en cuanto a los elementos n , tomamos su fila y la barajamos a través de los 4 elementos del registro a través de la instrucción de barajar , para multiplicar por la fila de la matriz m en el registro. Por ejemplo, para el elemento n._00 (recuerde que su desplazamiento en el registro tiene el índice 3), será:

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

En una forma simplificada, el algoritmo se ve así:

 //   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; 

Ahora cambiamos los elementos n en el algoritmo a la combinación aleatoria correspondiente, la multiplicación por _mm_mul_ps , la suma por _mm_add_ps , y ya está . Funciona Sin embargo, el código parece mucho peor de lo que parecía el algoritmo. Para este código, IACA emitió: x86 - 18.89, x64 - 16 ciclos . Esto es 4 veces más rápido que el anterior. En SSE registre el 4 ° flotador. Relación casi lineal.

Decorar implementación de SSE

Aún así, en el código se ve horrible. Intentaremos mejorar esto escribiendo un poco de azúcar sintáctica.


El compilador puede alinear perfectamente estas funciones (aunque a veces sin __forceinline de ninguna manera).


Y así, ya es mucho mejor y más legible. Para esto, IACA produjo aproximadamente el resultado esperado: x86 - 19 (y ¿por qué no fraccional?), X64 - 16 . De hecho, el rendimiento no ha cambiado, pero el código es mucho más hermoso y comprensible.

Poca contribución a la optimización futura

Presentemos una mejora más a nivel de una función que apareció recientemente en la versión de hierro. Operación de adición múltiple (fma) . fma (a, b, c) = a * b + c .


¿Por qué es esto necesario? En primer lugar, para una futura optimización. Por ejemplo, simplemente puede reemplazar mad con fma en el código terminado a través de las mismas macros que desee. Pero sentaremos las bases para la optimización ahora:


IACA: x86 - 18.89, x64 - 16 . De nuevo fraccional. Aún así, IACA a veces produce resultados extraños. El código no ha cambiado tanto. Probablemente incluso un poco peor. Pero la optimización a veces requiere tales sacrificios.

Pasamos al ahorro a través de _mm_stream

Varias guías de optimización recomiendan una vez más no extraer el caché para operaciones de guardado masivo. Esto generalmente se justifica cuando se procesan vértices que son miles o más. Pero para las matrices, esto quizás no sea tan importante. Sin embargo, lo agregaré de todos modos.


Nada ha cambiado en el tiempo aquí, de la palabra en absoluto. Pero, de acuerdo con las recomendaciones, ahora no volvemos a tocar el caché.

Implementación AVX

Opción base AVX

A continuación, pasamos a la siguiente etapa de optimización. El 4 ° flotador está incluido en el registro SSE, y en el AVX ya es 8. Es decir, existe una posibilidad teórica de reducir el número de operaciones realizadas y aumentar la productividad, si no a la mitad, y luego al menos 1,5 veces. Pero algo me dice que no todo será tan simple con la transición a AVX. ¿Podemos obtener los datos necesarios de registros dobles?

Tratemos de resolverlo. Nuevamente, escribimos nuestro algoritmo de multiplicación usado anteriormente. No podría hacer esto, pero es más conveniente lidiar con el código cuando todo está cerca y no tiene que desplazarse media página hacia arriba.

 //    : 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 

En la salida, esperamos obtener el resultado en ymm = {r0: r1} y ymm = {r2: r3} . Si en la versión SSE nuestro algoritmo se generalizó a columnas, ahora necesitamos generalizarlo a filas. Por lo tanto, actuar como en el caso de la opción SSE no funcionará.

Si consideramos la matriz m en los registros ymm , obtenemos ymm = {m0: m1} e ymm = {m2: m3}, respectivamente. Anteriormente, solo teníamos columnas matriciales en el registro, y ahora columnas y filas.

Si intenta actuar como antes, debe multiplicar ymm = {m0: m1} por el registro ymm = {n00, n00, n00, n00}: {n10, n10, n10, n10} . Dado que n00 y n01 están en la misma fila de la matriz n , a juzgar por el conjunto disponible de instrucciones AVX, dispersarlos por ymm será costoso. Tanto la combinación aleatoria como la permutación funcionan por separado para cada una de las cuatro patas de un flotador (alto y bajo xmm ) dentro de los registros ymm .

Si tomamos ymm de la matriz n , obtenemos ambos elementos n00 y n10 en el más alto de 2 xmm dentro del registro ymm . {n00, n10, n20, n30}: {n01, n11, n21, n31} . Por lo general, el índice para las instrucciones existentes es de 0 a 3. Y las direcciones flotan solo dentro de un registro xmm de dos dentro del registro ymm . No es posible transferir n10 del xmm más antiguo al más joven a bajo precio . Y luego este enfoque debe repetirse varias veces. No podemos soportar tal pérdida de medidas. Es necesario llegar a algo más.

Solíamos generalizar columnas, pero ahora filas. Por lo tanto, intentaremos ir de una manera un poco diferente. Necesitamos obtener el resultado en {r0: r1} . Esto significa que el algoritmo debe mejorarse no en líneas separadas del algoritmo, sino en dos a la vez. Y aquí, lo que fue una desventaja en el trabajo de barajar y permutar , será una ventaja para nosotros. Observamos lo que tendremos en los registros ymm cuando consideramos la matriz n .

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

Sí, notamos que en diferentes partes xmm del registro ymm tenemos los elementos 00 y 01 . Se pueden multiplicar mayúsculas y minúsculas mediante el comando permute en {_00, _00, _00, _00}: {_ 01, _01, _01, _01} , indicando solo un índice 3 para ambas partes xmm . Esto es exactamente lo que necesitamos. De hecho, los coeficientes también se usan en diferentes líneas. Solo ahora en el registro ymm correspondiente para la multiplicación será necesario mantener {m0: m0} , es decir, la primera fila duplicada de la matriz m .

Entonces, pintamos el algoritmo con más detalle. Leemos las filas dobles de la matriz m en registros ymm :

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

Y luego calcularemos la multiplicación 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) 

Reescribimos más 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> 

O en 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> 

Todo parece estar claro.


Aquí están los números interesantes de IACA: x86 - 12.53, x64 - 12 . Aunque, por supuesto, quería algo mejor. Perdí algo.

Optimización AVX más azúcar sintáctico

Parece que en el código anterior, AVX no se utilizó en todo su potencial. Encontramos que en lugar de establecer dos líneas idénticas en el registro ymm , podemos usar broadcast , que puede llenar el registro ymm con dos valores xmm idénticos. Además, en el camino, agregue algo de "azúcar sintáctico" para las funciones AVX.


Y aquí los resultados ya son más interesantes. IACA produce números: x86 - 10, x64 - 8.58 , que se ve mucho mejor, pero aún no 2 veces.

Opción AVX + FMA (final)

Hagamos un intento más. Ahora sería lógico recuperar el conjunto de instrucciones de FMA nuevamente, ya que se agregó a los procesadores después de AVX. Simplemente cambie mul + add individual para una operación. Aunque todavía utilizamos las instrucciones de multiplicación para dar al compilador más oportunidades de optimización, y al procesador para la ejecución paralela de multiplicaciones. Por lo general, miro el código generado en el ensamblador para asegurarme de qué opción es mejor.

En este caso, necesitamos calcular a * b + c * d + e * f + g * h . Puede hacer esta frente: fma (a, b, fma (c, d, fma (e, f, g * h))) . Pero, como vemos, es imposible realizar una operación aquí sin completar la anterior. Y esto significa que no podremos usar la capacidad de hacer multiplicaciones emparejadas, como la tubería SIMD nos permite hacer. Si transformamos ligeramente los cálculos fma (a, b, c * d) + fma (e, f, g * h) , veremos que podemos paralelizar los cálculos. Primero haga dos multiplicaciones independientes, y luego dos operaciones fma independientes.


IACA: x86 - 9.21, x64 - 8 . Ahora está muy bueno. Alguien probablemente dirá qué se puede hacer aún mejor, pero no sé cómo.

Puntos de referencia

Noto de inmediato que estas cifras no deben tomarse como la verdad última. Incluso con una prueba fija, nadan dentro de ciertos límites. Y aún más, se comportan de manera diferente en diferentes plataformas. Con cualquier optimización, tome medidas específicamente para su caso.

Tabla de contenidos

• Función: nombre de la función. La terminación en s funciona con una transmisión, mov diferente y normal (sin transmisión). Agregado para mayor claridad, ya que esto es lo suficientemente importante.
• Ciclos de IACA: número de ticks por función calculados por IACA
• Ciclos medidos: número medido de medidas (menos es más)
• Aceleración de IACA: número de medidas en una línea cero / número de medidas en una línea
• Aceleración medida: número de medidas en la línea cero / número de medidas en la línea (cuanto más, mejor)

Para loop_m, los ticks del artículo se multiplicaron por 64. Es decir, este es un valor muy aproximado. De hecho, resultó de esa manera.

i3-3770:

i7-8700K:

Código de prueba en fuente. Si hay sugerencias razonables sobre cómo mejorarlas, escriba los comentarios.

BONIFICACIÓN del reino de la ficción.

En realidad, es del reino de la ficción porque si vi procesadores con soporte para AVX512, tal vez en las imágenes. Sin embargo, intenté implementar el algoritmo. Aquí no explicaré nada, una analogía completa con AVX + FMA. El algoritmo es el mismo, solo que menos operaciones.


Los números son fantásticos: x86 - 4.79, x64 - 5.42 (IACA con arquitectura SKX). Esto a pesar del hecho de que el algoritmo tiene 64 multiplicaciones y 48 adiciones.

Código PS del artículo

Esta es mi primera experiencia escribiendo un artículo. Gracias a todos por sus comentarios. Ayudan a mejorar el código y el artículo.