为什么与 SSE 的矩阵乘法较慢？

Question

我有一个矩阵类（4x4）

class matrix {
public:
matrix() {}
matrix(float m11,float m21,float m31,float m41,
       float m12,float m22,float m32,float m42,
       float m13,float m23,float m33,float m43,
       float m14,float m24,float m34,float m44);
matrix(const float*);
matrix(const matrix&);

matrix operator *(const matrix& other)const;

static const matrix identity;
private:
union {
    float m[16];
    struct {
        float m11,m21,m31,m41;
        float m12,m22,m32,m42;
        float m13,m23,m33,m43;
        float m14,m24,m34,m44;
    };
    struct {
        float element[4][4];
    };
};
};

下面是乘法运算符的第一个实现，

matrix matrix::operator*(const matrix &other) const{
return matrix(
    m11*other.m11+m12*other.m21+m13*other.m31+m14*other.m41,
    m21*other.m11+m22*other.m21+m23*other.m31+m24*other.m41,
    m31*other.m11+m32*other.m21+m33*other.m31+m34*other.m41,
    m41*other.m11+m42*other.m21+m43*other.m31+m44*other.m41,
    m11*other.m12+m12*other.m22+m13*other.m32+m14*other.m42,
    m21*other.m12+m22*other.m22+m23*other.m32+m24*other.m42,
    m31*other.m12+m32*other.m22+m33*other.m32+m34*other.m42,
    m41*other.m12+m42*other.m22+m43*other.m32+m44*other.m42,
    m11*other.m13+m12*other.m23+m13*other.m33+m14*other.m43,
    m21*other.m13+m22*other.m23+m23*other.m33+m24*other.m43,
    m31*other.m13+m32*other.m23+m33*other.m33+m34*other.m43,
    m41*other.m13+m42*other.m23+m43*other.m33+m44*other.m43,
    m11*other.m14+m12*other.m24+m13*other.m34+m14*other.m44,
    m21*other.m14+m22*other.m24+m23*other.m34+m24*other.m44,
    m31*other.m14+m32*other.m24+m33*other.m34+m34*other.m44,
    m41*other.m14+m42*other.m24+m43*other.m34+m44*other.m44
);
}

我尝试使用 sse 指令来加速以下版本，

matrix matrix::operator*(const matrix &other) const{
float r[4][4];
__m128 c1=_mm_loadu_ps(&m11);
__m128 c2=_mm_loadu_ps(&m12);
__m128 c3=_mm_loadu_ps(&m13);
__m128 c4=_mm_loadu_ps(&m14);
for (int i = 0;i < 4; ++i) {
    __m128 v1 = _mm_set1_ps(other.element[i][0]);
    __m128 v2 = _mm_set1_ps(other.element[i][1]);
    __m128 v3 = _mm_set1_ps(other.element[i][2]);
    __m128 v4 = _mm_set1_ps(other.element[i][3]);

    __m128 col = _mm_add_ps(
        _mm_add_ps(_mm_mul_ps(v1,c1),_mm_mul_ps(v2,c2)),
        _mm_add_ps(_mm_mul_ps(v3,c3),_mm_mul_ps(v4,c4))
    );
    _mm_storeu_ps(r[i], col);
}
return matrix(&r[0][0]);
}

但在我的 macbookpro 上，第一个版本执行 100000 个矩阵乘法大约需要 6ms，第二个版本大约需要 8ms。 我想知道为什么会这样。 也许是因为 cpu 管道使第一个版本运行并发计算并且加载/保存滞后于第二个版本？

Answer 1

在第一种（标量）情况下，当您允许编译器按其认为最好的方式优化代码时，您将从大规模指令并行性中受益。通过arranging the code so as to minimize data dependencies，即使这可能导致需要更多的总指令，每条指令也可以在不同的执行单元上同时运行。有 很多 寄存器可用，因此大多数值可以保持注册状态，从而最大限度地减少昂贵的内存读取需求，即使在需要内存读取时，它们几乎可以免费完成，而其他操作正在完成，这要归功于乱序执行调度。我会进一步推测您正在从这里的 μ-op 缓存中受益，其好处是补偿了增加的代码大小。

在第二种（并行）情况下，您正在创建重要的数据依赖关系。即使编译器发出最佳目标代码（当您使用内在函数时不一定会出现这种情况），强制这种并行性也会产生成本。如果你ask the compiler to show you an assembly listing，你可以看到。在操作之间打包和重新排序 SSE 寄存器中的浮点操作数需要大量 shufps 指令。这在现代英特尔架构上只需要一个周期^*，但随后的addps 和mulps 操作无法并行执行。他们必须等待它完成。这段代码很有可能遇到了一个困难的 μ-op 吞吐量瓶颈。 （您可能还会在此代码中支付未对齐的数据损失，但这在现代架构中是最小的。）

换句话说，您已经用并行性（以更大的代码为代价）换取了增加的数据依赖性（尽管代码更小）。至少，这将是我的半受过教育的猜测，查看您的示例代码的反汇编。在这种情况下，您的基准非常清楚地告诉您它对您不利。

如果您指示编译器假定支持 AVX，情况可能会发生变化。如果目标体系结构不支持 AVX，编译器别无选择，只能将您的 _mm_set1_ps 内在转换为一对 movss、shufps 指令。如果启用 AVX 支持，您将获得一条 vbroadcastss 指令，这可能会更快，尤其是在支持 AVX2 的情况下，您可以从寄存器到寄存器进行广播（而不仅仅是从内存到寄存器）。借助 AVX 支持，您还可以获得 VEX 编码指令的好处。

---

_{^* 尽管在某些较旧的架构（如 Core 2）上，shufps 是基于整数的指令，因此当它后跟浮点指令时会导致延迟，例如addps 或 mulps。我不记得这个问题到底是什么时候修复的，但在 Sandy Bridge 及以后肯定不是问题。}