如何在 C++ 中编写矩阵乘法的快速代码？

Question

我为矩阵乘法编写了一些 C++ 代码。我使用vector<double> 来保存矩阵条目，并使用一系列3 个嵌套的for 循环来逐项计算乘法。事实证明，这非常慢（对于 900*500 和 500*500 的矩阵乘法，在我的 macbook air 上大约需要 10 秒）。是什么原因？是我使用了错误的矩阵表示还是代码中存在很大缺陷？

    for (int c_b=0;c_b<B.n_c;c_b++)
    {
        vector<double> vtmp(A.n_r);
        for (int r_a=0;r_a<A.n_r;r_a++)
        {
            sum=0;
            for (int i=0;i < A.n_c;i++)
            {
                sum=sum+A.mat[r_a+i*A.n_r]*B.mat[i+c_b*B.n_r];
            }
            vtmp[r_a]=sum;
        }
        Cvv[c_b]=vtmp;
    }

更新：此问题已通过在 Lapack 中使用子例程得到解决。

Answer 1

这里有一些关于性能改进的建议。

将向量移出循环
向量的创建需要时间。将声明移到任何 for 循环之前：

vector<double> vtmp(A.n_r);
for (int c_b=0;c_b<B.n_c;c_b++)
{
    for (int r_a=0;r_a<A.n_r;r_a++)
    {
        //...
    }
}

剖析分配计算。
如果没有任何分析或基准测试，赋值语句看起来占用的时间最多。将其分解为单独的步骤以帮助编译器，这样您就可以查看是否可以更优化地执行计算。
原文：

sum=sum+A.mat[r_a+i*A.n_r]*B.mat[i+c_b*B.n_r];

剖析[1]：

const int A_Index = r_a + i * A.n_r;
const int B_Index = i + c_b*B.n_r;
sum = sum + A.mat[A_Index] * B.mat[B_Index];

Dissected[2]（使用更多变量）：

const int temp1 = i * A.n_r;
const int temp2 = c_b * B.n_r;
const int A_Index = r_a + temp1;
const int B_Index = i + temp2;
sum = sum + A.mat[A_Index] * B.mat[B_Index];

以上内容可能有助于编译器选择最佳处理器指令。

使用局部变量
理想情况下，您希望处理器在重新加载之前在数据缓存中从矩阵中获取尽可能多的顺序位置。像这样的：

int ATemp1 = A[w];
int ATemp2 = A[x];
int ATemp3 = A[y];
int ATemp4 = A[z];

int BTemp1 = B[e];
int BTemp2 = B[f];
int BTemp3 = B[g];
int BTemp4 = B[h];

sum = sum + ATemp1 * BTemp1;
sum = sum + ATemp2 * BTemp2;
sum = sum + ATemp3 * BTemp3;
sum = sum + ATemp4 * BTemp4;