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?