c - 高效的 4x4 矩阵乘法(C 语言与汇编语言)

Question

我正在寻找一种用 C 语言将两个 4x4 矩阵相乘的更快、更复杂的方法。我目前的研究重点是具有 SIMD 扩展的 x86-64 汇编。到目前为止，我已经创建了一个比原始 C 实现快 6 倍的函数，这超出了我对性能改进的预期。不幸的是，只有在没有优化标志用于编译时(GCC 4.7)，这才会成立。使用 -O2，C 变得更快，我的努力变得毫无意义。

我知道现代编译器利用复杂的优化技术来实现近乎完美的代码，通常比巧妙的手工汇编更快。但在少数性能关键的情况下，人类可能会尝试与编译器争夺时钟周期。特别是，当可以探索一些以现代 ISA 为后盾的数学时(就像我的情况一样)。

我的函数如下所示(AT&T 语法，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

它通过处理封装在 128 位 SSE 寄存器中的四个 float ，在每次迭代中计算结果矩阵的整列。通过一些数学运算(操作重新排序和聚合)和 mullps/addps 4xfloat 包的并行乘法/加法指令，可以实现完全矢量化。代码重用用于传递参数的寄存器(%rdi、%rsi、%rdx:GNU/Linux ABI)，受益于(内部)循环展开并将一个矩阵完全保存在 XMM 寄存器中以减少内存读取。如您所见，我已经研究了这个主题，并花时间尽我所能地实现它。

征服我的代码的朴素 C 计算如下所示:

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

我研究了上述 C 代码的优化汇编输出，在 XMM 寄存器中存储 float 时，不涉及任何并行操作——仅涉及标量计算、指针算术和条件跳转。编译器的代码似乎不那么刻意，但它仍然比我预期的矢量化版本要快 4 倍左右更有效。我确信总体思路是正确的——程序员做类似的事情并获得有益的结果。但是这里有什么问题呢？是否有任何我不知道的寄存器分配或指令调度问题？你知道任何 x86-64 汇编工具或技巧来支持我与机器的战斗吗？

Answer 1

4x4矩阵乘法是64次乘法和48次加法。使用 SSE 这可以减少到 16 次乘法和 12 次加法(和 16 次广播)。以下代码将为您执行此操作。它只需要 SSE ( #include <xmmintrin.h> )。数组 A , B , 和 C需要 16 字节对齐。使用水平指令，例如 hadd (SSE3) 和 dpps (SSE4.1) 将是 less efficient (特别是 dpps )。我不知道循环展开是否有帮助。

void M4x4_SSE(float *A, float *B, float *C) {
    __m128 row1 = _mm_load_ps(&B[0]);
    __m128 row2 = _mm_load_ps(&B[4]);
    __m128 row3 = _mm_load_ps(&B[8]);
    __m128 row4 = _mm_load_ps(&B[12]);
    for(int i=0; i<4; i++) {
        __m128 brod1 = _mm_set1_ps(A[4*i + 0]);
        __m128 brod2 = _mm_set1_ps(A[4*i + 1]);
        __m128 brod3 = _mm_set1_ps(A[4*i + 2]);
        __m128 brod4 = _mm_set1_ps(A[4*i + 3]);
        __m128 row = _mm_add_ps(
                    _mm_add_ps(
                        _mm_mul_ps(brod1, row1),
                        _mm_mul_ps(brod2, row2)),
                    _mm_add_ps(
                        _mm_mul_ps(brod3, row3),
                        _mm_mul_ps(brod4, row4)));
        _mm_store_ps(&C[4*i], row);
    }
}