高效的 4x4 矩阵乘法（C 与汇编）

我正在寻找一种更快、更棘手的方法来用 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 寄存器中的四个浮点来计算每次迭代结果矩阵的整列。通过一些数学运算（操作重新排序和聚合）就可以实现完整的矢量化mullps/addps4xfloat 包的并行乘法/加法指令。该代码重用了用于传递参数的寄存器（%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 寄存器中存储浮点数时，不涉及任何并行操作– 只是标量计算、指针算术和条件跳转。编译器的代码似乎不那么刻意，但它仍然比我预计的矢量化版本快 4 倍左右，稍微有效一些。我确信总体想法是正确的——程序员做类似的事情并获得回报。但这里出了什么问题呢？是否存在我不知道的寄存器分配或指令调度问题？你知道任何 x86-64 组装工具或技巧来支持我与机器的战斗吗？

4x4 矩阵乘法是 64 次乘法和 48 次加法。使用 SSE，这可以减少到 16 次乘法和 12 次加法（以及 16 次广播）。以下代码将为您完成此操作。它只需要SSE（#include <xmmintrin.h>）。数组A, B, and C需要16字节对齐。使用水平指令，例如hadd（上证3）和dpps（SSE4.1）将是效率较低 https://stackoverflow.com/questions/14967969/efficient-4x4-matrix-vector-multiplication-with-sse-horizontal-add-and-dot-prod（尤其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);
    }
}