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

23
我正在寻找一种更快、更巧妙的方法来在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 汇编工具或技巧来支持我对抗机器吗?

1
这正是我所做的 - 我使用了一种替代计算方法来适应SSE问题。实际上,它是一个不同的算法。问题可能在于,现在我还必须在指令级别上进行优化,因为在专注于算法的同时,我可能已经引入了数据依赖问题,低效的内存访问模式或其他黑魔法。 - Krzysztof Abramowicz
1
如果您在 C 函数的指针参数中添加 restrict 限定符并使用 -O3 编译,GCC 将对其进行矢量化处理。如果没有 restrict 限定符,编译器必须假设输出矩阵可能与输入矩阵之一相同。 - caf
我还有一个使用SSE指令集的实现,即使使用-O2编译只能提高约10%的速度,而在禁用优化时几乎要慢两倍。由于它非常依赖于优化,所以我决定使用纯汇编来充分利用SSE方法。不幸的是,这导致了最慢的解决方案。也许我应该回到内部函数并在这个领域进行更多实验。 - Krzysztof Abramowicz
@KrzysztofAbramowicz,AVX已经发布了相当长的时间。您是否考虑使用AVX进行此操作?我可以添加一些代码,展示如何在4x4矩阵乘法中获得比SSE快两倍的速度。至少对于将4x4矩阵数组乘以固定矩阵来说是这样。 - Z boson
我在问题如何每个周期实现4个FLOPC代码循环性能中发现了类似的考虑。 - Krzysztof Abramowicz
显示剩余5条评论
5个回答

38
4x4矩阵乘法需要进行64次乘法和48次加法。使用SSE技术,这可以减少到16次乘法和12次加法(还有16次广播)。以下代码可以为您完成此操作。它仅需要SSE (#include <xmmintrin.h>)。数组ABC需要16字节对齐。使用水平指令,如hadd (SSE3)和dpps (SSE4.1)将会效率更低(特别是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);
    }
}

非常感谢您的回答。这段代码看起来比我之前使用 SSE intrinsics 进行矩阵乘法的实验要好。它还可以在 -O2 下生成更好的汇编代码,并且比我的代码运行速度稍快。但我仍然想知道为什么我不能用纯汇编至少达到相同的结果。 - Krzysztof Abramowicz
如果你正在使用GCC,为什么不使用“-O3”进行编译呢? - Z boson
也许是因为我一直被告知 -O3 引入了激进的优化技术,可能不会提高性能,但可能会引入额外的成本,例如在展开循环或内联函数时增加代码大小。但你是对的 - 首先是 -O3,然后是低级别的优化! :-) 幸运的是,在我的例子中,这并没有太大的区别。 - Krzysztof Abramowicz

14

有一种方法可以加速代码并超越编译器。 它不涉及任何复杂的流水线分析或深层次的代码微观优化(这并不意味着它不能进一步受益于这些)。 优化使用了三个简单的技巧:

  1. 函数现在是32字节对齐的(这显著提高了性能),

  2. 主循环以反向方式进行,这减少了与零测试(基于EFLAGS)的比较,

  3. 指令级地址算术证明比“外部”指针计算更快(即使它需要两倍的加法“在3/4情况下”)。 它缩短了循环体四条指令,并减少了其执行路径内的数据依赖关系。请参阅相关问题

此外,该代码使用相对跳转语法,可抑制符号重定义错误,该错误会在将其放置在asm语句中并使用-O3进行编译后尝试内联时发生()。

    .text
    .align 32                           # 1. function entry alignment
    .globl matrixMultiplyASM            #    (for a faster call)
    .type matrixMultiplyASM, @function
matrixMultiplyASM:
    movaps   (%rdi), %xmm0
    movaps 16(%rdi), %xmm1
    movaps 32(%rdi), %xmm2
    movaps 48(%rdi), %xmm3
    movq $48, %rcx                      # 2. loop reversal
1:                                      #    (for simpler exit condition)
    movss (%rsi, %rcx), %xmm4           # 3. extended address operands
    shufps $0, %xmm4, %xmm4             #    (faster than pointer calculation)
    mulps %xmm0, %xmm4
    movaps %xmm4, %xmm5
    movss 4(%rsi, %rcx), %xmm4
    shufps $0, %xmm4, %xmm4
    mulps %xmm1, %xmm4
    addps %xmm4, %xmm5
    movss 8(%rsi, %rcx), %xmm4
    shufps $0, %xmm4, %xmm4
    mulps %xmm2, %xmm4
    addps %xmm4, %xmm5
    movss 12(%rsi, %rcx), %xmm4
    shufps $0, %xmm4, %xmm4
    mulps %xmm3, %xmm4
    addps %xmm4, %xmm5
    movaps %xmm5, (%rdx, %rcx)
    subq $16, %rcx                      # one 'sub' (vs 'add' & 'cmp')
    jge 1b                              # SF=OF, idiom: jump if positive
    ret

这是我迄今为止见过的最快的x86-64实现。如果有更快的汇编代码,请提供答案并投票支持!


我在使用这个时遇到了麻烦。我使用以下签名从C中调用它: void abramowicz_MM4x4(float *A, float *B, float *C); 然后我在另一个文件中有相应的汇编代码以匹配gcc名称修饰: .globl _Z16abramowicz_MM4x4PfS_S _Z16abramowicz_MM4x4PfS_S: 调用返回了不正确的值。可能出了什么问题? - Praxeolitic
1
问题在于参数的顺序被颠倒了。对于任何想要尝试这个的人,可以在C函数签名中翻转A和B,或者在汇编语言中翻转rdi和rsi。 - Praxeolitic
有人有上面的英特尔汇编翻译吗? - Jason Nelson
我在这里写了一篇关于该主题的扩展博客文章:http://nfrechette.github.io/2017/04/13/modern_simd_matrix_multiplication/。我还将汇编版本翻译成了可由Visual Studio使用的内容,尽管对我的版本进行了一些微小的更改以保持其二进制精确性。不过,我的非汇编版本速度稍快! - Nicholas Frechette

4
Sandy Bridge及以上版本扩展了指令集以支持8个元素的向量算术。考虑这个实现。
struct MATRIX {
    union {
        float  f[4][4];
        __m128 m[4];
        __m256 n[2];
    };
};
MATRIX myMultiply(MATRIX M1, MATRIX M2) {
    // Perform a 4x4 matrix multiply by a 4x4 matrix 
    // Be sure to run in 64 bit mode and set right flags
    // Properties, C/C++, Enable Enhanced Instruction, /arch:AVX 
    // Having MATRIX on a 32 byte bundry does help performance
    MATRIX mResult;
    __m256 a0, a1, b0, b1;
    __m256 c0, c1, c2, c3, c4, c5, c6, c7;
    __m256 t0, t1, u0, u1;

    t0 = M1.n[0];                                                   // t0 = a00, a01, a02, a03, a10, a11, a12, a13
    t1 = M1.n[1];                                                   // t1 = a20, a21, a22, a23, a30, a31, a32, a33
    u0 = M2.n[0];                                                   // u0 = b00, b01, b02, b03, b10, b11, b12, b13
    u1 = M2.n[1];                                                   // u1 = b20, b21, b22, b23, b30, b31, b32, b33

    a0 = _mm256_shuffle_ps(t0, t0, _MM_SHUFFLE(0, 0, 0, 0));        // a0 = a00, a00, a00, a00, a10, a10, a10, a10
    a1 = _mm256_shuffle_ps(t1, t1, _MM_SHUFFLE(0, 0, 0, 0));        // a1 = a20, a20, a20, a20, a30, a30, a30, a30
    b0 = _mm256_permute2f128_ps(u0, u0, 0x00);                      // b0 = b00, b01, b02, b03, b00, b01, b02, b03  
    c0 = _mm256_mul_ps(a0, b0);                                     // c0 = a00*b00  a00*b01  a00*b02  a00*b03  a10*b00  a10*b01  a10*b02  a10*b03
    c1 = _mm256_mul_ps(a1, b0);                                     // c1 = a20*b00  a20*b01  a20*b02  a20*b03  a30*b00  a30*b01  a30*b02  a30*b03

    a0 = _mm256_shuffle_ps(t0, t0, _MM_SHUFFLE(1, 1, 1, 1));        // a0 = a01, a01, a01, a01, a11, a11, a11, a11
    a1 = _mm256_shuffle_ps(t1, t1, _MM_SHUFFLE(1, 1, 1, 1));        // a1 = a21, a21, a21, a21, a31, a31, a31, a31
    b0 = _mm256_permute2f128_ps(u0, u0, 0x11);                      // b0 = b10, b11, b12, b13, b10, b11, b12, b13
    c2 = _mm256_mul_ps(a0, b0);                                     // c2 = a01*b10  a01*b11  a01*b12  a01*b13  a11*b10  a11*b11  a11*b12  a11*b13
    c3 = _mm256_mul_ps(a1, b0);                                     // c3 = a21*b10  a21*b11  a21*b12  a21*b13  a31*b10  a31*b11  a31*b12  a31*b13

    a0 = _mm256_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 2, 2));        // a0 = a02, a02, a02, a02, a12, a12, a12, a12
    a1 = _mm256_shuffle_ps(t1, t1, _MM_SHUFFLE(2, 2, 2, 2));        // a1 = a22, a22, a22, a22, a32, a32, a32, a32
    b1 = _mm256_permute2f128_ps(u1, u1, 0x00);                      // b0 = b20, b21, b22, b23, b20, b21, b22, b23
    c4 = _mm256_mul_ps(a0, b1);                                     // c4 = a02*b20  a02*b21  a02*b22  a02*b23  a12*b20  a12*b21  a12*b22  a12*b23
    c5 = _mm256_mul_ps(a1, b1);                                     // c5 = a22*b20  a22*b21  a22*b22  a22*b23  a32*b20  a32*b21  a32*b22  a32*b23

    a0 = _mm256_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 3, 3, 3));        // a0 = a03, a03, a03, a03, a13, a13, a13, a13
    a1 = _mm256_shuffle_ps(t1, t1, _MM_SHUFFLE(3, 3, 3, 3));        // a1 = a23, a23, a23, a23, a33, a33, a33, a33
    b1 = _mm256_permute2f128_ps(u1, u1, 0x11);                      // b0 = b30, b31, b32, b33, b30, b31, b32, b33
    c6 = _mm256_mul_ps(a0, b1);                                     // c6 = a03*b30  a03*b31  a03*b32  a03*b33  a13*b30  a13*b31  a13*b32  a13*b33
    c7 = _mm256_mul_ps(a1, b1);                                     // c7 = a23*b30  a23*b31  a23*b32  a23*b33  a33*b30  a33*b31  a33*b32  a33*b33

    c0 = _mm256_add_ps(c0, c2);                                     // c0 = c0 + c2 (two terms, first two rows)
    c4 = _mm256_add_ps(c4, c6);                                     // c4 = c4 + c6 (the other two terms, first two rows)
    c1 = _mm256_add_ps(c1, c3);                                     // c1 = c1 + c3 (two terms, second two rows)
    c5 = _mm256_add_ps(c5, c7);                                     // c5 = c5 + c7 (the other two terms, second two rose)

                                                                    // Finally complete addition of all four terms and return the results
    mResult.n[0] = _mm256_add_ps(c0, c4);       // n0 = a00*b00+a01*b10+a02*b20+a03*b30  a00*b01+a01*b11+a02*b21+a03*b31  a00*b02+a01*b12+a02*b22+a03*b32  a00*b03+a01*b13+a02*b23+a03*b33
                                                //      a10*b00+a11*b10+a12*b20+a13*b30  a10*b01+a11*b11+a12*b21+a13*b31  a10*b02+a11*b12+a12*b22+a13*b32  a10*b03+a11*b13+a12*b23+a13*b33
    mResult.n[1] = _mm256_add_ps(c1, c5);       // n1 = a20*b00+a21*b10+a22*b20+a23*b30  a20*b01+a21*b11+a22*b21+a23*b31  a20*b02+a21*b12+a22*b22+a23*b32  a20*b03+a21*b13+a22*b23+a23*b33
                                                //      a30*b00+a31*b10+a32*b20+a33*b30  a30*b01+a31*b11+a32*b21+a33*b31  a30*b02+a31*b12+a32*b22+a33*b32  a30*b03+a31*b13+a32*b23+a33*b33
    return mResult;
}

.xmm[].ymm[] 可能是更好的联合成员名称。除此之外,看起来不错。虽然需要进行相当多的洗牌,但将其存储到内存中以便进行广播式加载可能会更好。(除非编译器将其"优化"回到洗牌操作...) - Peter Cordes
1
在Haswell及以后的处理器上,vbroadcastss ymm,[mem]是加载端口中的单个uop。在SnB/IvB上,它是一个加载+port5洗牌操作。但这仍然比vshufps + vperm2f128(或vinsertf128)的2个port5洗牌操作要快。 - Peter Cordes
1
哦,算了,你正在进行两个单独的车道广播,并且permute2f128在另一个操作数上。是的,看起来不错。使用-march=haswell,4个乘法/加法对折叠成FMA:https://godbolt.org/g/9uEbhR。嗯,那些`_mm256_permute2f128_ps(same,same, 0)是广播,但编译器没有将它们转换为vinsertf128`。这就是你可以通过Haswell的广播128负载来节省洗牌端口uops的地方。 - Peter Cordes

3

我想知道转置其中一个矩阵是否有益。

考虑如何计算以下两个矩阵的乘积...

A1 A2 A3 A4        W1 W2 W3 W4
B1 B2 B3 B4        X1 X2 X3 X4
C1 C2 C3 C4    *   Y1 Y2 Y3 Y4
D1 D2 D3 D4        Z1 Z2 Z3 Z4

这将导致...
dot(A,?1) dot(A,?2) dot(A,?3) dot(A,?4)
dot(B,?1) dot(B,?2) dot(B,?3) dot(B,?4)
dot(C,?1) dot(C,?2) dot(C,?3) dot(C,?4)
dot(D,?1) dot(D,?2) dot(D,?3) dot(D,?4)

将一行和一列做点积很麻烦。

如果我们在相乘之前转置第二个矩阵呢?

A1 A2 A3 A4        W1 X1 Y1 Z1
B1 B2 B3 B4        W2 X2 Y2 Z2
C1 C2 C3 C4    *   W3 X3 Y3 Z3
D1 D2 D3 D4        W4 X4 Y4 Z4

现在我们不再是对一行和一列进行点积,而是对两行进行点积。这可能更好地利用SIMD指令。
希望这有所帮助。

5
使用SSE时,几乎从不使用两个向量的点积。相反,您可以同时执行四个点积。与标量代码相同,但使用SIMD寄存器。例如,对于四分量向量,这意味着您需要执行4个 _mm_mul_ps 和 3个 _mm_add_ps,以得到四个点积结果。 - Z boson
1
@redrum:我明白了。到目前为止,我一直在使用“mulps”和“haddps”的组合来进行点积和矩阵、向量乘法。看来我还需要做更多的调整。 - Sparky
1
hadd 有时候会有用处,但在这种情况下不需要使用。我从来没有发现 dpps 有用过。 - Z boson
@user1095108,如果SIMD代码使用得当,它看起来几乎与标量情况相同。考虑3D点积w = x1x2 + y1y2 +z1z1。这里的变量可以是标量或SIMD向量,结果w可以是单个数字或一次多个点积的结果。每种情况使用相同数量的指令:3个乘法和2个加法。但标量情况使用标量指令,而SIMD情况使用SIMD指令。 - Z boson
1
@user1095108,我从未使用 _mm_dp_ps_mm_hadd_ps 来进行单个点积。我会尝试重新组织我的代码,以避免使用它们。请阅读这篇文章 http://www.cdl.uni-saarland.de/papers/leissa_vecimp_tr.pdf。但是英特尔一定有其创建 _mm_dp_ps 的原因。我之前看过他们的一份说明。如果您无法更改代码并且必须逐个计算一个点积,则 _mm_dp_ps 可能具有某些优势,但据我所知,这只是一种小的改进,远不能与同时计算四个点积所获得的4倍因子相提并论。您可以编写代码来测试这一点。 - Z boson
显示剩余4条评论

-3

显然,你可以一次从四个矩阵中获取项,并使用相同的算法同时乘以四个矩阵。


详细说明一下...它真的回答了问题吗? - peterh
我认为从四个不同的输入矩阵中收集元素,然后将其散布回四个不同的结果矩阵,速度不会比使用类似于OP自己答案中所做的load+broadcast更快。 - Peter Cordes

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接