Wydajne mnożenie macierzy 4x4 (C vs montaż)

Szukam szybszego i bardziej skomplikowanego sposobu pomnożenia dwóch macierzy 4x4 w C. Moje obecne badania skupiają się na montażu x86-64 z rozszerzeniami SIMD. Do tej pory stworzyłem funkcję, która jest około 6 razy szybsza niż naiwna implementacja C, która przekroczyła moje oczekiwania dotyczące poprawy wydajności. Niestety, pozostaje to prawdą tylko wtedy, gdy do kompilacji nie są używane flagi optymalizacji (GCC 4.7). Z-O2, C staje się szybszy i mój wysiłek staje się bez znaczenia.

Wiem, że nowoczesne kompilatory wykorzystują złożone techniki optymalizacji, aby uzyskać niemal idealny kod, zwykle szybszy niż pomysłowy kawałek ręcznie wykonanego montażu. Ale w mniejszej liczbie przypadków krytycznych dla wydajności, człowiek może próbować walczyć o cykle zegara z kompilatorem. Zwłaszcza, gdy niektóre matematyki poparte nowoczesnym ISA mogą być eksplorowane (tak jak w moim przypadku).

Moja funkcja wygląda następująco (składnia 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

Oblicza całą kolumnę wynikowej macierzy na iterację, przetwarzając cztery zmienne spakowane w 128-bitowych rejestrach SSE. Pełna wektoryzacja jest możliwa dzięki odrobinie matematyki (zmiana kolejności operacji i agregacja) imullps/addps instrukcje do równoległego mnożenia / dodawania pakietów 4xfloat. Kod ponownie wykorzystuje rejestry przeznaczone do przekazywania parametrów (%rdi, %rsi, %rdx : GNU / Linux ABI), korzysta z (wewnętrznej) pętli rozwijania i trzyma jedną macierz całkowicie w rejestrach XMM w celu zmniejszenia odczytów pamięci. Widzisz, zbadałem temat i nie spieszyłem się, aby wdrożyć go najlepiej, jak potrafię.

Obliczenia naiwne C podbijające mój kod wyglądają tak:

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

Zbadałem zoptymalizowany wynik montażu powyższego kodu C, który podczas przechowywania elementów pływających w rejestrach XMM,nie obejmuje żadnych równoległych operacji - tylko obliczenia skalarne, arytmetyczne wskaźniki i skoki warunkowe. Kod kompilatora wydaje się być mniej rozmyślny, ale nadal jest nieco bardziej efektywny niż moja wersja wektoryzacyjna, która ma być około 4 razy szybsza. Jestem pewien, że ogólna idea jest poprawna - programiści robią podobne rzeczy z satysfakcjonującymi wynikami. Ale co tu jest nie tak? Czy są jakieś problemy z rejestracją lub planowaniem instrukcji, o których nie wiem? Czy znasz jakieś narzędzia lub sztuczki montażowe x86-64, które wspierają moją walkę z maszyną?

questionAnswers(5)

yourAnswerToTheQuestion