使用矢量化 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

这一行：

Z = _mm_add_ps(M, N);

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

Answer 2

这有很多错误。

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

两个循环都前进了 4，因此内循环的主体一次处理旧标量循环的 16 个步骤。除了它没有，它确实 "four things".

而且它们不是正确的东西：

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 的列组合，所以BY 等 m1 行中的第二个条目乘以 m2 列中的第一个条目，这种情况不会发生。有很多不同的方式来排列那个计算，但是这里实现的方式不是重新排列，它计算了一些错误的产品并且跳过了很多必要的产品。

从 m2 加载一小行很方便，广播来自 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];

.. 但是把它放在循环中是不好的，只有当乘积的结果是应该在这些位置求和的数字时才有意义。这个 load/sum/store 东西在内循环中，因为内循环是 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 时不会。

使用矢量化 C++ 的矩阵乘法

Matrix Multiplication using vectorized c++

c++

sse

simd

matrix-multiplication