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

Question

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

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 [=11=]xB1, %%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 [=11=]x0A, %%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)。

Answer 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 内在函数。

Return 按值然后只是添加一个函数的问题：

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] = {}; 或类似的 - 这样您就可以默认获得 'zero-vector' 或默认构造函数：

constexpr AVector () = default;

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

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

Answer 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 [=10=]xB1, %%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 [=10=]x0A, %%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，因此它可能不会在片段运行之后将这些寄存器重新加载到它们以前的值。一些额外的 clobbers 也可能会丢失。

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

SSE inline assembly and possible g++ optimization bug

c++

sse

vector

matrix