矩阵向量乘法性能

计算科学 矩阵 数字
2021-12-24 15:56:24

我一直在了解缓存大小对代码性能的影响。我编写了一个小代码来查看在 MATLAB 中使用列主循环如何比使用行主循环更好，因为 MATLAB 将矩阵存储在 FORTRAN 等主列中。我还与 MATLAB 的内部乘法例程进行了比较。这是代码和结果：

% row major access
tic
for i=1:n
    for j=1:n
        b1(i)=b1(i)+A(i,j)*x(j);
    end
end
t1(count)=toc;

% column major access
tic
for j=1:n
    for i=1:n
        b2(i)=b2(i)+A(i,j)*x(j);
    end
end
t2(count)=toc;

% column major vector ops
tic
for i=1:n
    b3(i)=b3(i)+A(i,:)*x(:);
end
t3(count)=toc;

% row major vector ops
for j=1:n
    b4(:)=b4(:)+A(:,j)*x(j);
end
t4(count)=toc;

% MATLAB built in
tic
b5=A*x;
t5(count)=toc;

% double vop
tic
b6=A(:,:)*x(:);
t6(count)=toc;

在此处输入图像描述

列主循环比行主循环快，但是 MATLAB 怎么这么快呢？

3个回答

MATLAB 循环很慢。这是众所周知的。如果你使用 C++ 或 Julia 之类的东西，你会更接近（比如 5x-10x 左右的 IIRC）。所以这就是它超级疯狂的原因。

尽管如此，MATLAB 调用的算法实际上是BLAS库。具体来说，MATLAB 在后台调用了Intel MKL，因此它是一个极快的多线程代码。事实上，它可以做各种各样的事情。它首先将所有输入加载到连续的堆栈分配缓冲区中，然后它对指针进行一些操作以使事情移动得更快，SIMD 的一些事情，并且这一切都是通过并行性来完成的。尝试和理解是很疯狂的事情。如果您对OpenBLAS是开源的（而不是 MKL 不是）感到好奇，那么您可以看看并尝试了解发生了什么，但这太疯狂了！

编辑：

我最初说“然后通过模板进行操作，这样它实际上是一个命令而不是使用定义的天真矩阵乘法给出”。我想它实际上并没有在实践中使用 Strassen。此外，这是用于矩阵乘法，而不是像所问的问题那样专门用于矩阵向量乘法。对不起。O(nlog27)O(n3)

Matlab 可能正在移交给供应商 BLAS（MKL 等）。这些库与 Matlab 双 for 循环（线程并行、SIMD 内在函数、算法重新阻塞、工作）完全不同。

此外，无论如何，for 循环在 Matlab 中通常都很棘手。您可能会看到加速（在这两种情况下），只需以某种形式的矢量符号重写计算（通过编写A(:,j)/来处理整个列/行A(i,:)，那种东西）。

对其他两个（优秀）答案的补充：

理解 OpenBLAS 之类的源代码可能是一项艰巨的任务。作为替代起点，您可以通读

如何编写快速数字代码： Chellappa 等人的小介绍。

第 5 节处理 MMM（矩阵-矩阵乘法）运算。你从一个“朴素”的三循环 MMM 内核开始，检查它的内存数据访问模式。然后作者向您展示如何：

  1. 通过阻塞和缓冲优化缓存内存使用（第 5.1 节）
  2. 使用阻塞和循环展开在 CPU 和寄存器级别优化 MMM（第 5.2 节）
  3. 通过为您的 CPU 架构选择最佳参数（缓存/寄存器块大小等）来“调整”性能（第 5.3 节）。您可以使用ATLAS（第 5.4 节）之类的工具来执行此操作。

这不会让您接近高度优化的库（如英特尔 MKL）的性能，但您仍然应该看到相对于幼稚内核的显着改进。嘿，这是一个有趣的练习:)

其它你可能感兴趣的问题
上一篇有人用 Julia 写过 PDE 求解器吗？ 下一篇有限元代码的健全性检查雅可比