SIMD矩阵乘法_编程黑洞网

我最近开始玩弄SIMD，并想出了以下用于矩阵乘法的代码。

首先，我尝试使用SIMD来实现它，就像在SISD中一样，只是将SIMD用于点积之类的事情。对于每个特定的条目，它实际上要慢一些（仍然试图弄清楚这一点）。

经过一番思考，我意识到我可以对所得矩阵逐行进行计算，通过将这样的寄存器排成一行（每一行是一个SIMD寄存器，每列是x，y，z，w部分）：

带矩阵\ $ A \ $，\ $ B \ $和计算\ $ C = A * B \ $：

             A_00 * B_00                 A_00 * B_01
                 +                           +
            A_01 * B_10                 A_01 * B_11
                 +                           +
            A_02 * B_20                 A_02 * B_21
                 +                           +
            A_03 * B_30                 A_03 * B_31
                 =                           =
C_00 = Dot(A_Row0, B_Col0), C_01 = Dot(A_Row0, B_Col1), ...


            A_10 * B_00                 A_10 * B_01
                 +                           +
            A_11 * B_10                 A_11 * B_11
                 +                           +
                ...                         ...
C_10 = Dot(A_Row1, B_Col0), C_11 = Dot(A_Row1, B_Col1), ...

如果对这些事情有更多经验的人能告诉我，我离一个好的解决方案还有多远。

 __m128 BCx = _mm_load_ps((float*)&B.Row0);
__m128 BCy = _mm_load_ps((float*)&B.Row1);
__m128 BCz = _mm_load_ps((float*)&B.Row2);
__m128 BCw = _mm_load_ps((float*)&B.Row3);


// Calculate Row0 in resulting matrix
__m128 ARx = _mm_set1_ps(A.Row0.X);
__m128 ARy = _mm_set1_ps(A.Row0.Y);
__m128 ARz = _mm_set1_ps(A.Row0.Z);
__m128 ARw = _mm_set1_ps(A.Row0.W);

__m128 X = _mm_mul_ps(ARx, BCx);
__m128 Y = _mm_mul_ps(ARy, BCy);
__m128 Z = _mm_mul_ps(ARz, BCz);
__m128 W = _mm_mul_ps(ARw, BCw);

__m128 R = _mm_add_ps(X, _mm_add_ps(Y, _mm_add_ps(Z, W)));
_mm_storeu_ps((float*)&Result.Row0, R);

// Calculate Row1 in resulting matrix
ARx = _mm_set1_ps(A.Row1.X);
ARy = _mm_set1_ps(A.Row1.Y);
ARz = _mm_set1_ps(A.Row1.Z);
ARw = _mm_set1_ps(A.Row1.W);

X = _mm_mul_ps(ARx, BCx);
Y = _mm_mul_ps(ARy, BCy);
Z = _mm_mul_ps(ARz, BCz);
W = _mm_mul_ps(ARw, BCw);

R = _mm_add_ps(X, _mm_add_ps(Y, _mm_add_ps(Z, W)));
_mm_storeu_ps((float*)&Result.Row1, R);

// Calculate Row2 in resulting matrix
ARx = _mm_set1_ps(A.Row2.X);
ARy = _mm_set1_ps(A.Row2.Y);
ARz = _mm_set1_ps(A.Row2.Z);
ARw = _mm_set1_ps(A.Row2.W);

X = _mm_mul_ps(ARx, BCx);
Y = _mm_mul_ps(ARy, BCy);
Z = _mm_mul_ps(ARz, BCz);
W = _mm_mul_ps(ARw, BCw);

R = _mm_add_ps(X, _mm_add_ps(Y, _mm_add_ps(Z, W)));
_mm_storeu_ps((float*)&Result.Row2, R);

// Calculate Row3 in resulting matrix
ARx = _mm_set1_ps(A.Row3.X);
ARy = _mm_set1_ps(A.Row3.Y);
ARz = _mm_set1_ps(A.Row3.Z);
ARw = _mm_set1_ps(A.Row3.W);

X = _mm_mul_ps(ARx, BCx);
Y = _mm_mul_ps(ARy, BCy);
Z = _mm_mul_ps(ARz, BCz);
W = _mm_mul_ps(ARw, BCw);

R = _mm_add_ps(X, _mm_add_ps(Y, _mm_add_ps(Z, W)));
_mm_storeu_ps((float*)&Result.Row3, R);

使用完全优化（q431207 9q）在Visual Studio 2013编译器上，这大约是我典型的SISD版本的两倍（不确定要多少）。

这是我的SISD版本： br />

/Ox

SSE4上有一个点积说明可用：_mm_dp_ps。

嗨，谢谢您的回复，我也许应该提到我现在正以SSE2为目标，但我还没有完全研究目标受众可以合理预期的CPU类型。我也不确定sse4点产品的性能特征，有人提到它在某些硬件上速度较慢，如果有人能说明它会很棒，因为我的google-fu让我失败了。

@glampert：不过，DPPS效率不高。在Intel Haswell上，它具有14个周期的延迟，每个insn吞吐量有2个周期。 MULPS具有5个周期的延迟和每个insn吞吐量0.5个周期。（1/3的延迟，4倍的吞吐量）。参见agner.org/optimize。如果要执行大量操作，请垂直执行这些操作（fma，addps，mulps），最后只进行水平操作（haddps / dpps）。

这是C还是C ++？它们都是不同的语言。

它被编译为C ++，但更多地采用了C风格。问题主要是关于simd指令的正确用法和顺序，希望C和C ++之间的差异不会太大。随意假设有哪一个可以帮助您编写答案，并为显示两种语言的不同SIMD解决方案提供加分。

#1 楼

重组为函数以使其成为DRYer
现在，您的代码结构非常令人讨厌，而不是非常DRY，等等。我建议的第一件事是对其进行重组以使其在函数中起作用。我不确定您的Mat4结构是如何实现的，但是您在注释中指出它是连续的，因此我基于此假设。我建议将其封装到这种函数中。

void dotFourByFourMatrix(const Mat4* left, const Mat4* right, Mat4* result) {
    const __m128 BCx = _mm_load_ps((float*)&B.Row0);
    const __m128 BCy = _mm_load_ps((float*)&B.Row1);
    const __m128 BCz = _mm_load_ps((float*)&B.Row2);
    const __m128 BCw = _mm_load_ps((float*)&B.Row3);

    float* leftRowPointer = &left->Row0;
    float* resultRowPointer = &result->Row0;

    for (unsigned int i = 0; i < 4; ++i, leftRowPointer += 4, resultRowPointer += 4) {
        __m128 ARx = _mm_set1_ps(leftRowPointer[0]);
        __m128 ARy = _mm_set1_ps(leftRowPointer[1]);
        __m128 ARz = _mm_set1_ps(leftRowPointer[2]);
        __m128 ARw = _mm_set1_ps(leftRowPointer[3]);

        __m128 X = ARx * BCx;
        __m128 Y = ARy * BCy;
        __m128 Z = ARz * BCz;
        __m128 W = ARw * BCw;

        __m128 R = X + Y + Z + W;
        _mm_store_ps(resultRowPointer, R);
    }
}

您会注意到我对指针的增量进行了一些假设-我的假设是像Row0或X这样的成员变量仅仅是为了方便起见指向单个16成员连续数组的指针。如果那是不正确的，那将会中断。我也将类似BCx的变量标记为const，因为它们似乎好像不应该改变。可以在每个系统上完美运行-看起来leftRowPointer和resultRowPointer的增量是不确定的（在最后一次迭代之后），但是不太可能引起问题。我在StackOverfow上问了一个有关您可能有兴趣阅读的问题。
性能提升
SIMD通常超级容易做到，并且可以加快速度。只要您的数据布局良好（对齐，连续等），您就应该具有非常好的缓存行为，这是SIMD的最重要方面。真正改善此问题的唯一方法是添加预取功能（该文章很棒，并且在预取方面有很好的介绍），但是通常您必须稍微调整一下预取距离-即预取的距离。这是因为预取仍需要时间才能完成，因此，如果对下一个矩阵进行预取，除非当前计算花费足够的时间来掩盖预取，否则它可能不会提高任何速度。没有时间，就无法确定它是多少。预先知道您可以很好地处理哪种矩阵，预取不会给您任何好处。指针及其指向的任何内容）。再次，我没有看到任何证据，除非那是Row0之类的方法。 SIMD寄存器-使用_mm_set_ps通常表示您不知所措。我在这里看不到任何证据，但为将来参考可能有用。
使用更好和更少的内在函数
我将您的功能更改为不使用_mm_mul_ps或_mm_add_ps，而只使用运算符-更加易于阅读。坦率地说，我认为您可以在不牺牲太多可读性的情况下将更多内容压缩在一起-我不知道所有临时职位都是必需的。我还摆脱了对_mm_storeu_ps的调用，将其替换为对齐调用-这将大大加快速度，并且如果您有任何方法可以确保以理想的对齐方式分配矩阵。 >或将所有内在函数放在一起
总体而言，我发现内在函数有点难以阅读，而且它们的跨平台性也较低。我强烈建议使用Agner Fog的vectorclass库之类的库，以兼具可移植性和可读性。
如果您不能使用或不使用3rd party库，编写一个小包装类确实很容易，更具可读性，如果您喜欢模板和宏，也可以使其变得非常可移植。
演员
所有这些转换让我有些紧张-再次，我对它一无所知数据的格式，但似乎没有必要使用显式强制转换。如果必须强制转换，请不要使用C样式强制转换-我通常更喜欢static_cast。您可以在此处了解有关为何不应该在C ++代码中使用C样式强制转换的更多信息。总而言之，C样式强制转换将执行以下操作：

C强制转换是使用(type)object或type(object)进行强制转换。 C样式转换定义为以下第一个成功的转换：
const_cast static_cast（尽管忽略访问限制）static_cast（请参见上文），然后const_cast reinterpret_cast reinterpret_cast，然后const_cast

之所以很危险，是因为C样式的强制转换可能会变成reinterpret_cast（通过与上述相同的链接，增加了强调）：

reinterpret_cast是最危险的演员，应谨慎使用。它将一种类型直接转换为另一种类型-例如将值从一个指针转换为另一种指针，或将指针存储在int或其他各种讨厌的事物中。在很大程度上，使用reinterpret_cast获得的唯一保证是，通常，如果将结果强制转换回原始类型，则将获得完全相同的值（但如果中间类型小于原始类型，则不会得到相同的值）。

问题的大小和算法的选择
现在矩阵很小，如果保持不变，不要介意这一点，但是如果您认为它们会变大，则可能需要研究平铺乘法-与大型行x大行乘法相比，缓存行为要好得多。
按原样，使用这种大小的矩阵，除非您要进行大量处理，否则可能不会将会导致巨大的性能差异，而手动矢量化为您带来的可读性打击可能不那么值得。我从您的个人资料中看到您从事游戏开发，因此您可能做了很多，但是如果是这样，我建议您改为在GPU上进行。
DPPS
最后，我想突出显示注释中要使用_mm_dp_ps（DPPS）指令的一点。正如彼得·科德斯（Peter Cordes）指出的（也在评论中），这是一个坏主意。水平操作破坏了SIMD寄存器的全部内容-实际上，我前段时间实际上为点积的不同SIMD实现编写了一些时序测试，包括使用该指令，如果可以找到它们，我将上传结果图。如您所料，使用DPPS会大大降低您的速度。
编辑：您的评论指出这并不准确，但我将其保留在这里，因为这很重要。
我不知道您的矩阵是如何实现的，但是它们似乎不是连续的记忆，而是向量（或某些向量）的向量。这将损害您的缓存行为-尝试使用一维向量/数组，然后使用（假设以行为主）matrix[row*numColumns + col]进行访问。您应该看到缓存行为有所改善。

\ $ \ begingroup \ $
阅读代码会发现DPPS在我的代码中从未使用过，这在另一条评论中得到了建议，其他人则回答说它很慢。因为我只对SIMD代码感兴趣，并且可以修改其余的（但它们只是连续内存），所以未显示矩阵的实现（问题仅是SISD示例代码）。强制转换是由于将矩阵实现为包含4个结构（每行一个vector4）的结构。我将不得不对此投票，直到答案基于实际给出的代码。
\ $ \ endgroup \ $
–彼得
16年1月7日在21:20

\ $ \ begingroup \ $
@Peter我从来没有说过-我只是说清楚使用它是一个坏主意。请随意投票-我不打算在这个答案上做更多的事情。
\ $ \ endgroup \ $
–丹农
16年1月7日在21:22

\ $ \ begingroup \ $
“在评论中建议您使用DPPS-正如彼得·科德斯（Peter Cordes）指出的，这通常是一个坏主意。”建议我阅读注释，但不阅读显示的实际代码。
\ $ \ endgroup \ $
–彼得
16年1月7日在21:24

\ $ \ begingroup \ $
@Peter我希望其余答案表明我确实读过您的代码。
\ $ \ endgroup \ $
–丹农
16年1月7日在21:27

\ $ \ begingroup \ $
唯一的表明是您使用_mm_storeu_ps进行了评论，这是有效的观点。我仍然对您如何通过评论实际上不在代码中的内容来开始回答感到困惑。在这一点上，这个问题太老了，不再与我相关，我只是希望它对寻求良好实现的其他人有用。
\ $ \ endgroup \ $
–彼得
16年1月7日在21:30