[使用矢量化c ++的矩阵乘法

Question

我正在尝试编写C ++代码以使用SIMD进行矩阵乘法，但结果是错误的这是我的代码

    void mat_sse(DATA m1[][SIZE], DATA m2[][SIZE], DATA mout[][SIZE])
{

    DATA prod = 0;

    __m128 X, Y, Z, M, N;

    for(int i=0; i<SIZE; i=i+1){
    Z[0] = Z[1] = Z[2] = Z[3] = 0;
    for(int k=0; k< SIZE; k=k+4){

        for( int j=0; j<SIZE; j=j+4){
            X = _mm_load_ps(&m1[i][k]);
            Y = _mm_load_ps(&m2[k][j]);
            M = _mm_mul_ps(X, Y);
            Z = _mm_add_ps(M, N);
            mout[i][j] += Z[0];
        mout[i][j+1] += Z[1];
        mout[i][j+2] += Z[2];
        mout[i][j+3] += Z[3];
        }

    }

    }

    return ;

}

其中尺寸为const int SIZE = 40;你能帮忙吗？

Answer 1

这有很多问题。

for(int k=0; k< SIZE; k=k+4){
    for( int j=0; j<SIZE; j=j+4){

两个循环都前进4，因此内部循环的主体立即处理旧的标量循环的16步。除非没有，否则它会执行“四件事”。

它们不正确：

X = _mm_load_ps(&m1[i][k]);
Y = _mm_load_ps(&m2[k][j]);
M = _mm_mul_ps(X, Y);

因此，内部循环的每次迭代都从m1中提取相同的小行向量，并从m2中提取下一个小行向量，然后将它们逐点相乘。那不行例如，如果我们有两个4x4矩阵：（部分显示）

A B C D   X Y Z W
E . . .   S . . .
I . . . × T . . .
M . . .   U . . .

内部循环的迭代将计算AX，BY，CZ和DW。 AX确实应该出现在结果中，但真正的矩阵乘法不涉及BY：m1的行与m2的columns组合在一起，因此BY依此类推。 m1的行乘以m2列中的第一项，不会发生。有许多种方法可以安排该计算，但是此处实现的方法不是重新安排，它计算了一些错误的乘积，并且跳过了许多必要的乘积。

从m2加载一小行很方便，并且[[broadcast从m1加载单项很方便。这样，乘积在mout中仅占一行，因此可以对其进行累加并写入结果中，而无需进一步改组。

通过您已经完成的最后一部分的方式，

mout[i][j] += Z[0];
mout[i][j+1] += Z[1];
mout[i][j+2] += Z[2];
mout[i][j+3] += Z[3];

..但将其放在循环中是不好的，只有当乘积的结果是应该加到这些位置的数字时，才有意义。这种加载/求和/存储的事情在内部循环中，因为内部循环是j循环，但是可以通过交换j和k循环来解决：（未测试）
for (int i = 0; i < SIZE; i++) {
    for (int j = 0; j < SIZE; j += 4) {
        __m128 sum = _mm_setzero_ps();
        for (int k = 0; k < SIZE; k++) {
            __m128 entry = _mm_set1_ps(m1[i][k]);
            __m128 row  = _mm_load_ps(&m2[k][j]);
            sum = _mm_add_ps(sum, _mm_mul_ps(entry, row));
        }
        _mm_store_ps(&mout[i][j], sum);
    }
}

由于各种原因，该代码仍然很慢：

通过addps的循环传递依赖性比可用吞吐量慢。使用更多独立的累加器。
• 每个算术运算的负载过多。
• 对于中大型矩阵，请使用缓存阻止。不过size = 40时不行。

Answer 2

在此行：

Z = _mm_add_ps(M, N);

[N未初始化，因此Z将成为垃圾。