Eu estou procurando uma maneira mais rápida e mais complicada de multiplicar duas matrizes 4x4 em C. Minha pesquisa atual está focada na montagem x86-64 com extensões SIMD. Até agora, criei uma função que é cerca de 6x mais rápida que uma implementação ingênua de C, que excedeu minhas expectativas quanto à melhoria de desempenho. Infelizmente, isso permanece verdadeiro somente quando nenhum sinalizador de otimização é usado para compilação (GCC 4.7). Com-O2, C se torna mais rápido e meu esforço se torna sem sentido.

Eu sei que os compiladores modernos fazem uso de técnicas de otimização complexas para obter um código quase perfeito, geralmente mais rápido que uma peça engenhosa de montagem manual. Mas em uma minoria de casos críticos de desempenho, um humano pode tentar lutar por ciclos de clock com o compilador. Especialmente, quando alguma matemática apoiada com um ISA moderno pode ser explorada (como é no meu caso).

Minha função é a seguinte (sintaxe AT & T, GNU Assembler):

    .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

Ele calcula uma coluna inteira da matriz resultante por iteração, processando quatro flutuadores compactados em registradores SSE de 128 bits. A vectorização completa é possível com um pouco de matemática (reordenação e agregação de operação) emullps/addps instruções para multiplicação paralela / adição de pacotes 4xfloat. O código reutiliza registradores destinados a passar parâmetros (%rdi, %rsi, %rdx : ABI GNU / Linux), beneficia-se do desenrolamento de loop (interno) e mantém uma matriz inteiramente em registradores XMM para reduzir as leituras de memória. Você pode ver, pesquisei o assunto e aproveitei para implementá-lo da melhor maneira possível.

O cálculo C ingênuo conquistando meu código é assim:

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]);
}

Eu investiguei a saída de montagem otimizada do código C acima, que, enquanto armazena floats em registradores XMM,não envolve quaisquer operações paralelas - apenas cálculos escalares, aritmética de ponteiro e saltos condicionais. O código do compilador parece ser menos deliberado, mas ainda é um pouco mais efetivo do que minha versão vetorizada deve ser cerca de 4x mais rápida. Tenho certeza de que a idéia geral está correta - os programadores fazem coisas semelhantes com resultados recompensadores. Mas o que está errado aqui? Existe alguma alocação de registro ou problemas de programação de instrução dos quais não conheço? Você conhece alguma ferramenta ou truque de montagem x86-64 para apoiar minha batalha contra a máquina?