c++ - 为什么用 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 次矩阵乘法大约需要 6 毫秒，第二个版本大约需要 8 毫秒。我想知道为什么会这样。也许是因为 cpu 管道使第一个版本运行并发计算并且加载/保存滞后于第二个版本？

Answer 1

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

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

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

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

---

_{^* 尽管在某些较旧的架构(如 Core 2)上，shufps 是基于整数的指令，因此在其后跟一个指令时会导致延迟浮点指令，如 addps 或 mulps。我不记得这个问题是什么时候修复的，但肯定这在 Sandy Bridge 和之后的版本上不是问题。}