SSE 内联汇编和可能的 g++ 优化错误答案

【问题标题】：SSE inline assembly and possible g++ optimization bugSSE 内联汇编和可能的 g++ 优化错误
【发布时间】：2017-07-17 17:31:18
【问题描述】：

让我们从代码开始。我有两种结构，一种用于向量，另一种用于矩阵。

struct AVector
    {
    explicit AVector(float x=0.0f, float y=0.0f, float z=0.0f, float w=0.0f):
        x(x), y(y), z(z), w(w) {}
    AVector(const AVector& a):
        x(a.x), y(a.y), z(a.z), w(a.w) {}

    AVector& operator=(const AVector& a) {x=a.x; y=a.y; z=a.z; w=a.w; return *this;}

    float x, y, z, w;
    };

struct AMatrix
    {
    // Row-major
    explicit AMatrix(const AVector& a=AVector(), const AVector& b=AVector(), const AVector& c=AVector(), const AVector& d=AVector())
        {row[0]=a; row[1]=b; row[2]=c; row[3]=d;}
    AMatrix(const AMatrix& m) {row[0]=m.row[0]; row[1]=m.row[1]; row[2]=m.row[2]; row[3]=m.row[3];}

    AMatrix& operator=(const AMatrix& m) {row[0]=m.row[0]; row[1]=m.row[1]; row[2]=m.row[2]; row[3]=m.row[3]; return *this;}

    AVector row[4];
    };

接下来，对这些结构执行计算的代码。使用内联 ASM 和 SSE 指令的点积：

inline AVector AVectorDot(const AVector& a, const AVector& b)
    {
    // XXX
    /*const double v=a.x*b.x+a.y*b.y+a.z*b.z+a.w*b.w;

    return AVector(v, v, v, v);*/

    AVector c;

    asm volatile(
        "movups (%1), %%xmm0\n\t"
        "movups (%2), %%xmm1\n\t"
        "mulps %%xmm1, %%xmm0\n\t"          // xmm0 -> (a1+b1, , , )
        "movaps %%xmm0, %%xmm1\n\t"         // xmm1 = xmm0
        "shufps $0xB1, %%xmm1, %%xmm1\n\t"  // 0xB1 = 10110001
        "addps %%xmm1, %%xmm0\n\t"          // xmm1 -> (x, y, z, w)+(y, x, w, z)=(x+y, x+y, z+w, z+w)
        "movaps %%xmm0, %%xmm1\n\t"         // xmm1 = xmm0
        "shufps $0x0A, %%xmm1, %%xmm1\n\t"  // 0x0A = 00001010
        "addps %%xmm1, %%xmm0\n\t"          // xmm1 -> (x+y+z+w, , , )
        "movups %%xmm0, %0\n\t"
        : "=m"(c)
        : "r"(&a), "r"(&b)
        );

    return c;
    }

矩阵转置：

inline AMatrix AMatrixTranspose(const AMatrix& m)
    {
    AMatrix c(
        AVector(m.row[0].x, m.row[1].x, m.row[2].x, m.row[3].x),
        AVector(m.row[0].y, m.row[1].y, m.row[2].y, m.row[3].y),
        AVector(m.row[0].z, m.row[1].z, m.row[2].z, m.row[3].z),
        AVector(m.row[0].w, m.row[1].w, m.row[2].w, m.row[3].w));

    // XXX
    /*printf("AMcrix c:\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n",
        c.row[0].x, c.row[0].y, c.row[0].z, c.row[0].w,
        c.row[1].x, c.row[1].y, c.row[1].z, c.row[1].w,
        c.row[2].x, c.row[2].y, c.row[2].z, c.row[2].w,
        c.row[3].x, c.row[3].y, c.row[3].z, c.row[3].w);*/

    return c;
    }

矩阵-矩阵乘法 - 转置第一个矩阵，因为当我将它存储为列主要，第二个作为行主要时，我可以使用点积执行乘法。

inline AMatrix AMatrixMultiply(const AMatrix& a, const AMatrix& b)
    {
    AMatrix c;

    const AMatrix at=AMatrixTranspose(a);

    // XXX
    /*printf("AMatrix at:\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n",
        at.row[0].x, at.row[0].y, at.row[0].z, at.row[0].w,
        at.row[1].x, at.row[1].y, at.row[1].z, at.row[1].w,
        at.row[2].x, at.row[2].y, at.row[2].z, at.row[2].w,
        at.row[3].x, at.row[3].y, at.row[3].z, at.row[3].w);*/

    for(int i=0; i<4; ++i)
        {
        c.row[i].x=AVectorDot(at.row[0], b.row[i]).w;
        c.row[i].y=AVectorDot(at.row[1], b.row[i]).w;
        c.row[i].z=AVectorDot(at.row[2], b.row[i]).w;
        c.row[i].w=AVectorDot(at.row[3], b.row[i]).w;
        }

    return c;
    }

现在是主要（双关语）部分的时间：

int main(int argc, char *argv[])
    {
    AMatrix a(
        AVector(0, 1, 0, 0),
        AVector(1, 0, 0, 0),
        AVector(0, 0, 0, 1),
        AVector(0, 0, 1, 0)
        );

    AMatrix b(
        AVector(1, 0, 0, 0),
        AVector(0, 2, 0, 0),
        AVector(0, 0, 3, 0),
        AVector(0, 0, 0, 4)
        );

    AMatrix c=AMatrixMultiply(a, b);

    printf("AMatrix c:\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n    [%5.2f %5.2f %5.2f %5.2f]\n",
        c.row[0].x, c.row[0].y, c.row[0].z, c.row[0].w,
        c.row[1].x, c.row[1].y, c.row[1].z, c.row[1].w,
        c.row[2].x, c.row[2].y, c.row[2].z, c.row[2].w,
        c.row[3].x, c.row[3].y, c.row[3].z, c.row[3].w);

    AVector v(1, 2, 3, 4);
    AVector w(1, 1, 1, 1);

    printf("Dot product: %f (1+2+3+4 = 10)\n", AVectorDot(v, w).w);

    return 0;
    }

在上面的代码中，我创建了两个矩阵，将它们相乘并打印出结果矩阵。如果我不使用任何编译器优化（g++ main.cpp -O0 -msse），它工作正常。启用优化 (g++ main.cpp -O1 -msse) 结果矩阵为空（所有字段为零）。取消注释任何标有 XXX 的块会使程序写入正确的结果。

在我看来，GCC 优化了 AMatrixMultiply 函数中的矩阵，因为它错误地认为它没有在使用 SSE 内联编写的 AVectorDot 中使用。

最后几行检查点积函数是否真的有效，是的，确实有效。

所以，问题是：我做错了什么或理解错了什么，或者这是 GCC 中的某种错误？我的猜测是以上 7:3 的混合。

我使用的是 GCC 版本 5.1.0 (tdm-1)。

【问题讨论】：

这对我来说在clang中有效。
只使用内在函数。不仅编译器对它们的语义有更好的理解（理论上），你的__asm__ 也有错误。例如，它不知道 xmm 寄存器正在被破坏。你可能会侥幸逃脱，这在某些方面更糟。例如，如果你真的想做水平操作（不适合 SIMD），可以使用 _mm_dp_ps，甚至可以使用宏 _MM_TRANSPOSE4_PS。也可以在您的数据成员前面尝试alignas (16)。

标签： c++ matrix vector sse

【解决方案1】：

这也是使用 SSE 乘以矩阵的一种非常低效的方法。如果它比现代 CPU 上具有如此多浮点吞吐量的标量实现要快得多，我会感到惊讶。这里概述了一种更好的方法，不需要显式转置：

AMatrix & operator *= (AMatrix & m0, const AMatrix & m1)
{
    __m128 r0 = _mm_load_ps(& m1[0][x]);
    __m128 r1 = _mm_load_ps(& m1[1][x]);
    __m128 r2 = _mm_load_ps(& m1[2][x]);
    __m128 r3 = _mm_load_ps(& m1[3][x]);

    for (int i = 0; i < 4; i++)
    {
        __m128 ti = _mm_load_ps(& m0[i][x]), t0, t1, t2, t3;

        t0 = _mm_shuffle_ps(ti, ti, _MM_SHUFFLE(0, 0, 0, 0));
        t1 = _mm_shuffle_ps(ti, ti, _MM_SHUFFLE(1, 1, 1, 1));
        t2 = _mm_shuffle_ps(ti, ti, _MM_SHUFFLE(2, 2, 2, 2));
        t3 = _mm_shuffle_ps(ti, ti, _MM_SHUFFLE(3, 3, 3, 3));

        ti = t0 * r0 + t1 * r1 + t2 * r2 + t3 * r3;
        _mm_store_ps(& m0[i][x], ti);
    }

    return m0;
}

在 gcc 和 clang 等现代编译器上，t0 * r0 + t1 * r1 + t2 * r2 + t3 * r3 实际上是在 __m128 类型上运行的；不过如果你愿意，你可以用 _mm_mul_ps 和 _mm_add_ps 内部函数替换它们。

按值返回只需添加如下函数：

inline AMatrix operator * (const AMatrix & m0, const AMatrix & m1)
{
    AMatrix lhs (m0); return (lhs *= m1);
}

就我个人而言，我只是将 float x, y, z, w; 替换为 alignas (16) float _s[4] = {}; 或类似名称 - 这样您就可以默认获得“零向量”或默认构造函数：

constexpr AVector () = default;

以及不错的构造函数，例如：

constexpr Vector (float x, float y, float z, float w)
        : _s {x, y, z, w} {}

【讨论】：

好吧，我想它效率低下，因为它是在我使用 SSE 和内联 ASM 时编写的。当我使用内在函数时，我在内存对齐方面遇到了一些问题，但我会再试一次。至于结构，我会坚持我的 x、y、z、w 字段——它们比构造函数更常用；
@crueltear - gcc 手册不如教程好。我找到的最好的内联汇编是here。在我给出的示例中，您始终可以使用 loadu 和 storeu 等效项。

【解决方案2】：

您的内联程序集缺少一些约束：

asm volatile(
    "movups (%1), %%xmm0\n\t"
    "movups (%2), %%xmm1\n\t"
    "mulps %%xmm1, %%xmm0\n\t"          // xmm0 -> (a1+b1, , , )
    "movaps %%xmm0, %%xmm1\n\t"         // xmm1 = xmm0
    "shufps $0xB1, %%xmm1, %%xmm1\n\t"  // 0xB1 = 10110001
    "addps %%xmm1, %%xmm0\n\t"          // xmm1 -> (x, y, z, w)+(y, x, w, z)=(x+y, x+y, z+w, z+w)
    "movaps %%xmm0, %%xmm1\n\t"         // xmm1 = xmm0
    "shufps $0x0A, %%xmm1, %%xmm1\n\t"  // 0x0A = 00001010
    "addps %%xmm1, %%xmm0\n\t"          // xmm1 -> (x+y+z+w, , , )
    "movups %%xmm0, %0\n\t"
    : "=m"(c)
    : "r"(&a), "r"(&b)
    );

GCC 不知道这个汇编程序片段会破坏%xmm0 和%xmm1，因此它可能不会在片段运行后将这些寄存器重新加载到它们之前的值。一些额外的破坏者也可能会丢失。

【讨论】：

嗯，是的，这可能是我正在经历的一部分。猜猜是时候告别手写 SSE 并回到我更了解的东西了。