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

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 的 block 会使程序写入正确的结果。

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

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

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

我使用的是 GCC 版本 5.1.0 (tdm-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] = {}; 或类似的- 所以你默认得到一个“零 vector ”，或者一个默认的构造函数:

constexpr AVector () = default;

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

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

关于c++ - SSE 内联汇编和可能的 g++ 优化错误，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/45150313/