Ich suche nach einer schnelleren und schwierigeren Methode, um zwei 4x4-Matrizen in C zu multiplizieren. Meine aktuelle Forschung konzentriert sich auf die x86-64-Assemblierung mit SIMD-Erweiterungen. Bisher habe ich eine Funktion erstellt, die ungefähr 6x schneller ist als eine naive C-Implementierung, was meine Erwartungen an die Leistungsverbesserung übertroffen hat. Leider bleibt dies nur dann erhalten, wenn für die Kompilierung keine Optimierungsflags verwendet werden (GCC 4.7). Mit-O2, C wird schneller und meine Anstrengung wird bedeutungslos.

Ich weiß, dass moderne Compiler komplexe Optimierungstechniken einsetzen, um einen nahezu perfekten Code zu erzielen, der in der Regel schneller ist als ein geniales Stück handgefertigter Baugruppe. In einigen wenigen leistungskritischen Fällen kann ein Mensch jedoch versuchen, mit dem Compiler um Taktzyklen zu kämpfen. Insbesondere, wenn einige mit einer modernen ISA unterstützte Mathematik untersucht werden können (wie in meinem Fall).

Meine Funktion sieht wie folgt aus (AT & T-Syntax, 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

Es berechnet eine ganze Spalte der resultierenden Matrix pro Iteration, indem vier Floats verarbeitet werden, die in 128-Bit-SSE-Registern gepackt sind. Die vollständige Vektorisierung ist mit ein wenig Mathe (Umordnen und Aggregieren von Operationen) und möglichmullps/addps Anleitung zur parallelen Multiplikation / Addition von 4xfloat-Paketen. Der Code verwendet Register zur Übergabe von Parametern (%rdi, %rsi, %rdx : GNU / Linux ABI), profitiert vom (inneren) Loop-Unrolling und hält eine Matrix vollständig in XMM-Registern, um Speicherlesevorgänge zu reduzieren. Wie Sie sehen, habe ich das Thema recherchiert und mir die Zeit genommen, es so gut wie möglich umzusetzen.

Die naive C-Berechnung, die meinen Code erobert, sieht folgendermaßen aus:

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

Ich habe die optimierte Assembly-Ausgabe des obigen C-Codes untersucht, der beim Speichern von Floats in XMM-Registernbeinhaltet keine parallelen Operationen - Nur Skalarberechnungen, Zeigerarithmetik und bedingte Sprünge. Der Code des Compilers scheint weniger gewollt zu sein, aber er ist immer noch etwas effektiver als meine vektorisierte Version, von der erwartet wird, dass sie etwa 4x schneller ist. Ich bin sicher, dass die allgemeine Idee richtig ist - Programmierer machen ähnliche Dinge mit lohnenden Ergebnissen. Aber was stimmt hier nicht? Gibt es Probleme mit der Registerzuordnung oder der Befehlsplanung, die mir nicht bekannt sind? Kennen Sie irgendwelche x86-64-Montagewerkzeuge oder -Tricks, um meinen Kampf gegen die Maschine zu unterstützen?