【发布时间】：2013-09-01 05:44:20
【问题描述】：
我正在寻找一种在 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/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 组装工具或技巧来支持我与机器的战斗吗？
【问题讨论】：
-
最近的编译器可以比人类更好地进行微优化。专注于算法优化！
-
这正是我所做的——我使用了另一种计算来适应 SSE 的问题。它实际上是一种不同的算法。问题可能是，现在我还必须在指令级别对其进行优化，因为在专注于算法的同时，我可能引入了数据依赖性问题、无效的内存访问模式或其他一些黑魔法。
-
您最好使用通过
<immintrin.h>获得的 SSE 内在函数 - 您可以尝试其他方法，例如_mm_dp_ps和_MM_TRANSPOSE4_PS，而无需维护程序集。 -
如果将
restrict限定符添加到C 函数的指针参数并使用-O3进行编译，GCC 会将其向量化。如果没有restrict限定符，编译器必须假定输出矩阵可能与输入矩阵之一相同。 -
@BrettHale，我同意内在函数是这样做的方法，但
_mm_dp_ps或_MM_TRANSPOSE4_PS效率低下。查看我的回答和stackoverflow.com/questions/14967969/…
标签： c optimization assembly sse matrix-multiplication