用于矩阵计算的操作系统之间的巨大性能差异

Question

在我的两台计算机上，我尝试了以下代码：

N <- 10e3
M <- 2000
X <- matrix(rnorm(N * M), N)
system.time(crossprod(X))

第一个是标准笔记本电脑，此操作需要 1.7 秒。

> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

第二台是一台相当不错的台式电脑，耗时 17 秒。

> sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18.3

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

台式电脑的性能比笔记本电脑高，但这种矩阵计算需要多 10 倍的时间。

问题是否来自使用的默认 BLAS/LAPACK？

Answer 1

tldr： CentOS 使用单线程 OpenBLAS，Linux Mint 默认使用 Reference BLAS，但可以使用其他 BLAS 版本。

EPEL 提供的 CentOS R 软件包依赖于openblas-Rblas。这似乎是为 R 提供 BLAS 的 OpenBLAS 构建。因此，虽然看起来使用了 R 的 BLAS，但它实际上是 OpenBLAS。 LAPACK 版本始终是 R 提供的版本。

在 Debian 和 Mint 等衍生发行版上，r-base-core 取决于

• libblas3 | libblas.so.3
• liblapack3 | liblapack.so.3

默认情况下，这些由参考实现 libblas3 和 liblapack3 提供。这些并不是特别快，但您可以通过安装像 libopenblas-base 这样的软件包轻松替换它们。您可以通过update-alternatives 控制系统上使用的 BLAS 和 LAPACK。

为了控制 OpenBLAS 的线程数，我通常使用RhpcBLASctl：

N <- 20000
M <- 2000
X <- matrix(rnorm(N * M), N)
RhpcBLASctl::blas_set_num_threads(2)
system.time(crossprod(X))
#>        User      System verstrichen 
#>       2.492       0.331       1.339
RhpcBLASctl::blas_set_num_threads(1)
system.time(crossprod(X))
#>        User      System verstrichen 
#>       2.319       0.052       2.316

由于某种原因，从 R 中设置 environment variables OPENBLAS_NUM_THREADS、GOTO_NUM_THREADS 或 OMP_NUM_THREADS 并没有达到预期的效果。在 CentOS 上，即使 RhpcBLASctl 也无济于事，因为使用的 OpenBLAS 是单线程的。

Answer 2

R 随默认 BLAS 实现一起分发，但可能未针对您的计算机进行优化。 在 R as advised on the Installation guide of R 之前通过 ATLAS 或 OpenBLAS 安装 BLAS 的优化版本是要走的路。如果您单击Download R for Linux，然后单击debian/。据说：

您可能需要安装自动调整的 Atlas 或多线程 OpenBlas 库，以获得更高的线性代数运算性能

R 的源代码可以在here 下载，BLAS 实现位于R-3.5.0/src/extra/blas。例如，矩阵-矩阵乘法 dgemm 的 Fortran 源代码位于 blas.f 中，与大多数 BLAS 例程一起（在单个文件中！）。

函数的cmets指定：

 -- Written on 8-February-1989.
 Jack Dongarra, Argonne National Laboratory.
 Iain Duff, AERE Harwell.
 Jeremy Du Croz, Numerical Algorithms Group Ltd.
 Sven Hammarling, Numerical Algorithms Group Ltd.

可以在例程 dgemm 的 netlib 实现中找到完全相同的行。

相反，OpenBLAS 提供了不同的实现，每种处理器都有一个。参见例如this file devoted to dgemm for haswell microarchitecture。调用prefetcht0 进行预取，调用vfmadd231pd，这是一个矢量FMA SIMD 指令，一次执行4 次双精度d=a*b+c。

使用优化的 BLAS 可以节省时间。例如参见this benchmark，其中 netlib 的 dgemm() 持续 64 秒，而 MKL、OpenBLAS 或 ATLAS dgemm 只需不到 4 秒。

R 内部 BLAS 的情况可能比经典的 Netlib 库更糟糕。确实，如附录 R Installation and Administration 中的 A.3.1.5 Shared BLAS 所述：

R 提供了将 BLAS 编译为存储在 R_HOME/lib 中的动态库 libRblas 的选项，并将 R 本身和所有附加包链接到该库。 .... 使用共享 BLAS 可能存在性能劣势。 ...但是，实验表明，在许多情况下，只要使用高水平的编译器优化，使用共享 BLAS 的速度一样快。

看R的config.site文件，写的是g77/gfortran的优化级别是'-O2'。因此，如果 fortran 编译器不是 g77/gfortran，则调整 FFLAGS 选项可能很有用。在配置步骤中，应该有类似checking whether we are using the GNU Fortran 77 compiler... yes 的行（配置文件的第7521 行）。