Multiplicación eficiente de matrices 4x4 (C vs ensamblaje)

Multiplicación eficiente de matrices 4x4 (C vs ensamblaje)


Estoy buscando una forma más rápida y complicada de multiplicar dos matrices 4x4 en C. Mi investigación actual se centra en el ensamblaje x86-64 con extensiones SIMD. Hasta ahora, he creado una función que es aproximadamente 6 veces más rápida que una implementación de C ingenua, lo que ha superado mis expectativas de mejora del rendimiento. Desafortunadamente, esto se mantiene solo cuando no se usan banderas de optimización para la compilación (GCC 4.7). Con -O2 , C se vuelve más rápido y mi esfuerzo pierde sentido.


Sé que los compiladores modernos utilizan técnicas de optimización complejas para lograr un código casi perfecto, generalmente más rápido que una ingeniosa pieza de ensamblaje hecha a mano. Pero en una minoría de casos críticos para el rendimiento, un ser humano puede intentar luchar por los ciclos de reloj con el compilador. Especialmente, cuando se pueden explorar algunas matemáticas respaldadas con un ISA moderno (como es en mi caso).


Mi función tiene el siguiente aspecto (sintaxis de AT&T, ensamblador GNU):


    .text
.globl matrixMultiplyASM
.type matrixMultiplyASM, @function
matrixMultiplyASM:
movaps (%rdi), %xmm0 # fetch the first matrix (use four registers)
movaps 16(%rdi), %xmm1
movaps 32(%rdi), %xmm2
movaps 48(%rdi), %xmm3
xorq %rcx, %rcx # reset (forward) loop iterator
.ROW:
movss (%rsi), %xmm4 # Compute four values (one row) in parallel:
shufps $0x0, %xmm4, %xmm4 # 4x 4FP mul's, 3x 4FP add's 6x mov's per row,
mulps %xmm0, %xmm4 # expressed in four sequences of 5 instructions,
movaps %xmm4, %xmm5 # executed 4 times for 1 matrix multiplication.
addq $0x4, %rsi
movss (%rsi), %xmm4 # movss + shufps comprise _mm_set1_ps intrinsic
shufps $0x0, %xmm4, %xmm4 #
mulps %xmm1, %xmm4
addps %xmm4, %xmm5
addq $0x4, %rsi # manual pointer arithmetic simplifies addressing
movss (%rsi), %xmm4
shufps $0x0, %xmm4, %xmm4
mulps %xmm2, %xmm4 # actual computation happens here
addps %xmm4, %xmm5 #
addq $0x4, %rsi
movss (%rsi), %xmm4 # one mulps operand fetched per sequence
shufps $0x0, %xmm4, %xmm4 # |
mulps %xmm3, %xmm4 # the other is already waiting in %xmm[0-3]
addps %xmm4, %xmm5
addq $0x4, %rsi # 5 preceding comments stride among the 4 blocks
movaps %xmm5, (%rdx,%rcx) # store the resulting row, actually, a column
addq $0x10, %rcx # (matrices are stored in column-major order)
cmpq $0x40, %rcx
jne .ROW
ret
.size matrixMultiplyASM, .-matrixMultiplyASM

Calcula una columna completa de la matriz resultante por iteración, procesando cuatro flotantes empaquetados en registros SSE de 128 bits. La vectorización completa es posible con un poco de matemática (operación de reordenación y agregación) y mullps /addps instrucciones para la multiplicación/suma paralela de paquetes 4xfloat. El código reutiliza registros destinados a pasar parámetros (%rdi , %rsi , %rdx :GNU/Linux ABI), se beneficia del desenrollado de bucles (internos) y mantiene una matriz completamente en registros XMM para reducir las lecturas de memoria. Como puede ver, investigué el tema y me tomé mi tiempo para implementarlo lo mejor que pude.


El cálculo C ingenuo que conquista mi código se ve así:


void matrixMultiplyNormal(mat4_t *mat_a, mat4_t *mat_b, mat4_t *mat_r) {
for (unsigned int i = 0; i < 16; i += 4)
for (unsigned int j = 0; j < 4; ++j)
mat_r->m[i + j] = (mat_b->m[i + 0] * mat_a->m[j + 0])
+ (mat_b->m[i + 1] * mat_a->m[j + 4])
+ (mat_b->m[i + 2] * mat_a->m[j + 8])
+ (mat_b->m[i + 3] * mat_a->m[j + 12]);
}

He investigado la salida de ensamblado optimizada del código C anterior que, mientras almacena flotantes en registros XMM, no implica ninguna operación paralela – solo cálculos escalares, aritmética de punteros y saltos condicionales. El código del compilador parece ser menos deliberado, pero sigue siendo un poco más efectivo que mi versión vectorizada, que se esperaba que fuera aproximadamente 4 veces más rápida. Estoy seguro de que la idea general es correcta:los programadores hacen cosas similares con resultados gratificantes. Pero, ¿qué está mal aquí? ¿Hay algún problema de asignación de registros o programación de instrucciones que desconozco? ¿Conoces alguna herramienta o truco de ensamblaje x86-64 para apoyar mi batalla contra la máquina?


Respuestas:


Hay una forma de acelerar el código y superar al compilador. No implica ningún análisis de canalización sofisticado ni microoptimización de código profundo (lo que no significa que no pueda beneficiarse más de ellos). La optimización utiliza tres trucos simples:



  1. La función ahora está alineada en 32 bytes (lo que mejoró significativamente el rendimiento),


  2. El bucle principal va a la inversa, lo que reduce la comparación con una prueba cero (basado en EFLAGS),


  3. La aritmética de direcciones a nivel de instrucción demostró ser más rápida que el cálculo de punteros "externos" (aunque requiere el doble de adiciones «en 3/4 casos»). Acortó el cuerpo del bucle en cuatro instrucciones y redujo las dependencias de datos dentro de su ruta de ejecución. Ver pregunta relacionada.



Además, el código usa una sintaxis de salto relativa que suprime el error de redefinición de símbolo, que ocurre cuando GCC intenta alinearlo (después de colocarlo dentro de asm declaración y compilado con -O3 ).


    .text
.align 32 # 1. function entry alignment
.globl matrixMultiplyASM # (for a faster call)
.type matrixMultiplyASM, @function
matrixMultiplyASM:
movaps (%rdi), %xmm0
movaps 16(%rdi), %xmm1
movaps 32(%rdi), %xmm2
movaps 48(%rdi), %xmm3
movq $48, %rcx # 2. loop reversal
1: # (for simpler exit condition)
movss (%rsi, %rcx), %xmm4 # 3. extended address operands
shufps $0, %xmm4, %xmm4 # (faster than pointer calculation)
mulps %xmm0, %xmm4
movaps %xmm4, %xmm5
movss 4(%rsi, %rcx), %xmm4
shufps $0, %xmm4, %xmm4
mulps %xmm1, %xmm4
addps %xmm4, %xmm5
movss 8(%rsi, %rcx), %xmm4
shufps $0, %xmm4, %xmm4
mulps %xmm2, %xmm4
addps %xmm4, %xmm5
movss 12(%rsi, %rcx), %xmm4
shufps $0, %xmm4, %xmm4
mulps %xmm3, %xmm4
addps %xmm4, %xmm5
movaps %xmm5, (%rdx, %rcx)
subq $16, %rcx # one 'sub' (vs 'add' & 'cmp')
jge 1b # SF=OF, idiom: jump if positive
ret

Esta es la implementación x86-64 más rápida que he visto hasta ahora. ¡Apreciaré, votaré y aceptaré cualquier respuesta que proporcione una pieza de ensamblaje más rápida para ese propósito!