Eficiente multiplicación 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 ingenua de C, que ha superado mis expectativas para la mejora del rendimiento. Desafortunadamente, esto se mantiene así solo cuando no se utilizan indicadores 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 hacen uso de complejas técnicas de optimización para lograr un código casi perfecto, generalmente más rápido que una ingeniosa pieza de ensamblaje hecho a mano. Pero en una minoría de casos de rendimiento crítico, un humano puede intentar luchar por ciclos de reloj con el compilador. Especialmente, cuando se pueden explorar algunas matemáticas respaldadas con una ISA moderna (como en mi caso).
Mi función tiene el aspecto siguiente (sintaxis de AT&T, ensamblador de 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 flotadores empaquetados en registros SSE de 128 bits. La vectorización completa es posible con un poco de matemática (reordenación de la operación y agregación) ymullps
/addps
Instrucciones para la multiplicación / adición paralela de paquetes 4xfloat. El código reutiliza registros destinados a pasar parámetros (%rdi
, %rsi
, %rdx
: GNU / Linux ABI), se beneficia del desenrollado del bucle (interno) y mantiene una matriz completamente en los registros XMM para reducir las lecturas de memoria. Como puede ver, he investigado el tema y me tomé mi tiempo para implementarlo lo mejor que pueda.
El ingenuo cálculo de C 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 ensamblaje optimizada del código C anterior que, al almacenar los flotadores en los registros XMM,No implica operaciones paralelas. - Solo cálculos escalares, aritmética de punteros y saltos condicionales. El código del compilador parece ser menos deliberado, pero aún es un poco más efectivo de lo que se espera que mi versión vectorizada sea 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 registro o de programación de instrucciones que no conozca? ¿Conoces alguna herramienta de montaje x86-64 o trucos para apoyar mi batalla contra la máquina?