使用 C 和 OpenMP 进行矩阵乘法

Question

我正在尝试创建一种有效的算法，该算法可以将具有双精度的大值矩阵相乘。我已经创建了算法并首先在小矩阵上对其进行了测试；尝试后即 A{4096x4096}, B{4096x4096} 循环需要永远结束；例如，这两个矩阵生成 AB 花了我的电脑 30 多分钟才能完成。

我的电脑不是老旧的……它是六核 i7，我想对于桌面工作站来说还不错。在尺寸高达 1024x1024 的小矩阵上，它的完成速度相对较快，不到 30-40 秒，对于 2048x2048 大约需要 5 分钟……对于 16384x16384，它没有在 15 分钟内完成，我停止了执行……

我做错了什么还是可以预料到的？ :)

提前致谢！

代码如下：

/* calculate */
for(travx = 0; travx < m; travx++) {
    for(travy = 0; travy < n; travy++) {
        /* we only need to calculate it ourside of Z loop */
        tIndex = (travy)+(travx*n); 
        for(travz = 0; travz < p; travz++)
            {
                if(n==1)
                    {bIndex = ((n-1)*travy)+travz;
                     aIndex = ((p)*travx)+travz;} 
                else
                    {bIndex = ((n)*travz)+travy;
                     aIndex = ((p)*travx)+travz;}

                temp = atab_ptr[aIndex]*btab_ptr[bIndex];
                outtab_ptr[tIndex] =  outtab_ptr[tIndex] + temp;
            }
    }
}

这真的很简单......并且在小矩阵上给出了很好的结果......不知道如何在 10 秒内将双打相乘，尤其是在 p4 上......听起来有点可疑......特别是如果你考虑到 O (3) 问题的复杂性。

更新...根据反馈我调整了代码并且...嗯，主要是我对其进行了简化，小矩阵完成得更快，即 1024x1024 在 3 秒内完成，但 4096x4096 在 6 秒内完成分钟...修改后的代码是这样的：

for(travx = 0; travx < m; travx++) {
    for(travy = 0; travy < n; travy++) {
      for(travz = 0; travz < p; travz++)
        {outtab_ptr[travy+travx*n] = outtab_ptr[travy+travx*n] + atab_ptr[travy+p*travz] *  btab_ptr[travz+travx*p];}
    }
  }

Answer 1

如果可以的话，BLAS 是最好的选择。

话虽如此，从根本上说，矩阵乘法受到复杂性的限制，因此您必须更加聪明才能大幅缩短时间。矩阵是否以任何方式结构化？他们是tridiagonal 还是banded？他们是triangular 还是symmetric？

Answer 2

您的“高效”算法实际上效率很低。看看n 不是 1 时会发生什么：

bIndex = ((n)*travz)+travy;
aIndex = ((p)*travx)+travz;
temp = atab_ptr[aIndex]*btab_ptr[bIndex];

最里面的循环结束了travz，所以aIndex随着travz的每一个增量增加一步1。另一方面，bIndex 以n 的步长增加。因此，您正在访问 btab_ptr 的元素，这些元素在内存中不相邻，因此不在同一缓存行中。

更不用说最内层循环中的条件对可能的向量化有什么影响。

因此，如果所有矩阵的数据都可以放入 Core i7 的 L3 缓存中，那么您的算法运行速度可以接受，但一旦不是这种情况，您的性能就会急剧下降。然后进一步乘以 O(N^3) 复杂度。

Answer 3

好吧，矩阵乘法的简单方法是 O(n^3)。这意味着将两个矩阵相乘所需的时间随着输入的大小以三次方式增长。还有更有效的方法。在这里你可以看看：

http://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Matrix_algebra

这些方法仍然没有低于 O(n^2)。因此，随着矩阵大小的增加，完成时间以超线性方式越来越多，这是正常的。

话虽如此，您观察的时间是否过多，这取决于许多因素（您的机器、您的代码等）。

顺便说一句，您可以查看thread，其中提出了一个非常相似的问题。而且，除非您是出于教育目的，否则最好使用优化的库，例如 ATLAS。

在这里，您还有一个经典的document，关于如何优化应用程序以更好地使用内存。在该文档中，作者使用了对齐和预取等多种技术来优化矩阵乘法的性能。