我正在寻找一种更快、更复杂的方法来在 C 中将两个 4x4 矩阵相乘。我目前的研究主要集中在 x86-64 与 SIMD 扩展的汇编上。到目前为止，我已经创建了一个函数，它比简单的 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/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 寄存器中时，不涉及任何并行操作——仅涉及标量计算、指针算术和条件跳转。编译器的代码似乎没有那么刻意，但它仍然比我的矢量化版本更有效，预计速度大约快 4 倍。我确信一般的想法是正确的——程序员做类似的事情并获得有益的结果。但这里有什么问题？是否有任何我不知道的寄存器分配或指令调度问题？你知道任何 x86-64 组装工具或技巧来支持我与机器的战斗吗？