多核机器上单精度数组与双精度数组的矩阵乘法性能下降

Question

更新

不幸的是，由于我的疏忽，我有一个旧版本的 MKL (11.1) 链接到 numpy。较新版本的 MKL (11.3.1) 在 C 和从 python 调用时提供相同的性能。

令人费解的是，即使将已编译的共享库与较新的 MKL 显式链接，并通过 LD_* 变量指向它们，然后在 python 中执行 import numpy，以某种方式使 python 调用旧的 MKL 库。只有将 python lib 文件夹中的所有 libmkl_*.so 替换为更新的 MKL，我才能在 python 和 C 调用中匹配性能。

背景/图书馆信息。

矩阵乘法是通过 sgemm（单精度）和 dgemm（双精度）Intel 的 MKL 库调用，通过 numpy.dot 函数完成的。库函数的实际调用可以通过例如验证。 oprof。

这里使用 2x18 核心 CPU E5-2699 v3，因此总共有 36 个物理核心。 KMP_AFFINITY=分散。运行于 linux。

TL;DR

1) 为什么 numpy.dot，即使它调用相同的 MKL 库函数，与 C 编译代码相比最多慢两倍？

2) 为什么通过 numpy.dot 你会随着内核数量的增加而降低性能，而在 C 代码中却没有观察到相同的效果（调用相同的库函数）。

问题

我观察到在 numpy.dot 中进行 single/double 精度浮点数的矩阵乘法，以及直接从已编译的 C 共享库中调用 cblas_sgemm/dgemm 与从纯 C 代码内部调用相同的 MKL cblas_sgemm/dgemm 函数相比，性能明显更差。

import numpy as np
import mkl
n = 10000
A = np.random.randn(n,n).astype('float32')
B = np.random.randn(n,n).astype('float32')
C = np.zeros((n,n)).astype('float32')

mkl.set_num_threads(3); %time np.dot(A, B, out=C)
11.5 seconds
mkl.set_num_threads(6); %time np.dot(A, B, out=C)
6 seconds
mkl.set_num_threads(12); %time np.dot(A, B, out=C)
3 seconds
mkl.set_num_threads(18); %time np.dot(A, B, out=C)
2.4 seconds
mkl.set_num_threads(24); %time np.dot(A, B, out=C)
3.6 seconds
mkl.set_num_threads(30); %time np.dot(A, B, out=C)
5 seconds
mkl.set_num_threads(36); %time np.dot(A, B, out=C)
5.5 seconds

与上面完全相同，但使用双精度 A、B 和 C，您将得到： 3核：20s，6核：10s，12核：5s，18核：4.3s，24核：3s，30核：2.8s，36核：2.8s.

单精度浮点速度的提升似乎与缓存未命中有关。对于 28 核运行，这里是 perf 的输出。对于单精度：

perf stat -e task-clock,cycles,instructions,cache-references,cache-misses ./ptestf.py
631,301,854 cache-misses # 31.478 % of all cache refs

和双精度：

93,087,703 cache-misses # 5.164 % of all cache refs

C 共享库，用

编译

/opt/intel/bin/icc -o comp_sgemm_mkl.so -openmp -mkl sgem_lib.c -lm -lirc -O3 -fPIC -shared -std=c99 -vec-report1 -xhost -I/opt/intel/composer/mkl/include

#include <stdio.h>
#include <stdlib.h>
#include "mkl.h"

void comp_sgemm_mkl(int m, int n, int k, float *A, float *B, float *C);

void comp_sgemm_mkl(int m, int n, int k, float *A, float *B, float *C)
{
    int i, j;
    float alpha, beta;
    alpha = 1.0; beta = 0.0;

    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                m, n, k, alpha, A, k, B, n, beta, C, n);
}

Python包装函数，调用上面的编译库：

def comp_sgemm_mkl(A, B, out=None):
    lib = CDLL(omplib)
    lib.cblas_sgemm_mkl.argtypes = [c_int, c_int, c_int, 
                                 np.ctypeslib.ndpointer(dtype=np.float32, ndim=2), 
                                 np.ctypeslib.ndpointer(dtype=np.float32, ndim=2),
                                 np.ctypeslib.ndpointer(dtype=np.float32, ndim=2)]
    lib.comp_sgemm_mkl.restype = c_void_p
    m = A.shape[0]
    n = B.shape[0]
    k = B.shape[1]
    if np.isfortran(A):
        raise ValueError('Fortran array')
    if m != n:
        raise ValueError('Wrong matrix dimensions')
    if out is None:
        out = np.empty((m,k), np.float32)
    lib.comp_sgemm_mkl(m, n, k, A, B, out)

但是，通过 C 语言编译的二进制文件显式调用 MKL 的 cblas_sgemm / cblas_dgemm，并通过 C 中的 malloc 分配数组，与 python 相比，性能提高了近 2 倍代码，即 numpy.dot 调用。此外，未观察到性能随内核数量增加而下降的影响。 单精度矩阵乘法的最佳性能是 900 毫秒，这是在通过 mkl_set_num_cores 和运行将 C 代码与 numactl 结合使用所有 36 个物理内核时实现的 - -交错=全部。

对于 profiling/inspecting/understanding 这种情况，也许还有什么奇特的工具或建议？任何阅读 material 也非常感谢。

更新按照@Hristo Iliev 的建议，运行ning numactl --interleave=all ./ipython 没有改变时间（在噪音范围内），但提高了纯 C 二进制运行次。

Answer 1

我怀疑这是由于不幸的线程调度造成的。我能够重现与您类似的效果。 Python 在 ~2.2 秒时运行宁，而 C 版本显示出 1.4-2.2 秒的巨大变化。

正在申请： KMP_AFFINITY=scatter,granularity=thread 这可确保 28 个线程始终运行在同一个处理器线程上运行。

将运行时间减少到更稳定，C 为 ~1.24 秒，python 为 ~1.26 秒。

这是在 28 核双插槽 Xeon E5-2680 v3 系统上。

有趣的是，在非常相似的 24 核双插槽 Haswell 系统上，python 和 C 的性能几乎相同，即使没有线程关联/固定。

为什么python会影响排程？好吧，我假设它周围有更多运行时间环境。底线是，如果不固定，您的性能结果将是不确定的。

您还需要考虑，英特尔 OpenMP 运行时间会产生一个额外的管理线程，这可能会混淆调度程序。固定有更多选择，例如 KMP_AFFINITY=compact - 但由于某种原因，这在我的系统上完全混乱了。您可以将 ,verbose 添加到变量以查看运行time 如何固定您的线程。

likwid-pin 是一个有用的替代方案，提供更方便的控制。

一般来说，单精度应该至少和双精度一样快。双精度可能会更慢，因为：

双精度需要更多 memory/cache 带宽。
您可以构建单精度吞吐量更高的 ALU，但这通常不适用于 CPU，而是适用于 GPU。

我认为一旦你摆脱了性能异常，这就会反映在你的数字中。

当您增加 MKL/*gemm 的线程数时，请考虑

内存/共享缓存带宽可能成为瓶颈，限制可扩展性
Turbo 模式会在提高利用率时有效降低核心频率。即使您运行处于标称频率，这也适用：在 Haswell-EP 处理器上，AVX 指令将施加较低的 "AVX base frequency" - 但当使用较少的内核/热余量可用时，处理器可以超过该频率通常在短时间内甚至更多。如果您想要完全中性的结果，则必须使用 AVX 基本频率，对您来说是 1.9 GHz。已记录 here, and explained in one picture.

我不认为有一种真正简单的方法可以衡量您的应用程序如何受到不良计划的影响。您可以使用 perf trace -e sched:sched_switch 来公开它，并且有 some software 来可视化它，但这会带来很高的学习曲线。然后再次 - 对于并行性能分析，您应该固定线程。

多核机器上单精度数组与双精度数组的矩阵乘法性能下降

Performance degradation of matrix multiplication of single vs double precision arrays on multi-core machine

c

python

numpy

openmp

intel-mkl