时间戳

首页 > 解决方案 > 确定两台机器之间的 ddot 差异

scipy - 确定两台机器之间的 ddot 差异

问题描述

np.dot我目前有两台机器，它们为两个向量上的实例产生不同的输出。在没有深入挖掘从 NumPy 到 BLAS 的许多抽象层的情况下，我能够重现scipy.linalg.blas.ddot. 具体来说，考虑以下示例：

import numpy as np
from scipy.linalg.blas import ddot

u = np.array([0.13463703107579461093, -0.07773272613450200874, -0.98784132994666418170])
v = np.array([-0.86246572448831815283, -0.03715105562531360872, -0.50475010960748223354])
a = np.dot(v, u)
b = v[0]*u[0] + v[1]*u[1] + v[2]*u[2]
c = ddot(v, u)
print(f'{a:.25f}')
print(f'{b:.25f}')
print(f'{c:.25f}')

这会产生以下输出：

                      Machine 1                   Machine 2
a   0.3853810478481685120044631 0.3853810478481685675156143
b   0.3853810478481685120044631 0.3853810478481685120044631
c   0.3853810478481685120044631 0.3853810478481685675156143

同样，下面的 Cython 也产生了同样的差异：

cimport scipy.linalg.cython_blas
cimport numpy as np
import numpy as np

cdef np.float64_t run_test(np.double_t[:] a, np.double_t[:] b):
    cdef int ix, iy, n
    ix = iy = 1
    n = 3
    return scipy.linalg.cython_blas.ddot(&n, &a[0], &ix, &b[0], &iy)

a = np.array([0.13463703107579461093, -0.07773272613450200874, -0.98784132994666418170])
b = np.array([-0.86246572448831815283, -0.03715105562531360872, -0.50475010960748223354])
print(f'{run_test(a, b):.25f}')

所以，我试图了解是什么导致了这种情况。

有问题的机器分别运行 Windows 10（Intel(R) Core(TM) i7-5600U）和 Windows Server 2016（Intel(R) Xeon(R) Gold 6140）。

在这两种情况下，我都设置了新的 conda 环境，只有numpy, scipy,cython和它们的依赖项。我已经在环境上运行校验和，以确保最终包含的二进制文件同意并验证输出np.__config__.show()匹配。同样，我检查mkl.get_version_string()了两台机器上的输出是否一致。

这使我认为问题可能出在硬件差异上。我没有研究最终执行了哪些指令（缺乏在 Windows/MSVC 上调试 Cython 代码的直接方法），但我检查了两台机器是否都支持 AVX2/FMA，这似乎是差异的一个来源。

另一方面，我确实发现这两台机器支持不同的指令集。具体来说

          Machine 1 (i7)       Machine 2 (Xeon)
AVX                    Y                      Y
AVX2                   Y                      Y
AVX512CD               N                      Y
AVX512ER               N                      N
AVX512F                N                      Y
AVX512PF               N                      N
FMA                    Y                      Y

但是，我不知道确定这本身是否足以解释差异的好方法，或者它是否是一个红鲱鱼（？）

所以我的问题变成了：

从上面开始，有哪些自然的步骤可以尝试确定差异的原因？是组装时间，还是有更明显的东西？

标签： scipyblasintel-mklavx2avx512

解决方案

鉴于对该问题的出色评论，很明显支持的指令集之间的差异最终是罪魁祸首，实际上我们可以在运行 Cython 脚本时使用ListDLLs来发现 MKL 基于这两种情况加载不同的库。

对于 i7（机器 1）：

>listdlls64 python.exe | wsl grep mkl
0x00000000b9ff0000  0xe7e000   [...]\miniconda3\envs\npfloattest\Library\bin\mkl_rt.dll
0x00000000b80e0000  0x1f05000  [...]\miniconda3\envs\npfloattest\Library\bin\mkl_intel_thread.dll
0x00000000b3b40000  0x43ba000  [...]\miniconda3\envs\npfloattest\Library\bin\mkl_core.dll
0x00000000b0e50000  0x2ce5000  [...]\miniconda3\envs\npfloattest\Library\bin\mkl_avx2.dll
0x00000000b01f0000  0xc58000   [...]\miniconda3\envs\npfloattest\Library\bin\mkl_vml_avx2.dll
0x00000000f88c0000  0x7000     [...]\miniconda3\envs\npfloattest\lib\site-packages\mkl\_mklinit.cp37-win_amd64.pyd
0x00000000afce0000  0x22000    [...]\miniconda3\envs\npfloattest\lib\site-packages\mkl\_py_mkl_service.cp37-win_amd64.pyd

对于至强（机器 2）：

0x0000000057ec0000  0xe7e000   [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_rt.dll
0x0000000055fb0000  0x1f05000  [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_intel_thread.dll
0x0000000051bf0000  0x43ba000  [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_core.dll
0x000000004e1a0000  0x3a4a000  [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_avx512.dll
0x000000005c6c0000  0xc03000   [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_vml_avx512.dll
0x0000000079a70000  0x7000     [...]\Miniconda3\envs\npfloattest\lib\site-packages\mkl\_mklinit.cp37-win_amd64.pyd
0x000000005e830000  0x22000    [...]\Miniconda3\envs\npfloattest\lib\site-packages\mkl\_py_mkl_service.cp37-win_amd64.pyd

这非常强烈地表明，对 AVX512CD/AVX512F 的支持足以促使 MKL 使用不同的库，从而可能最终使用不同的指令集。

现在看看这实际上是如何展开的很有趣：发出了哪些指令，以及在具体数字示例中这意味着什么。

首先，让我们编写等效的 VC++ 程序来了解最终运行的指令：

typedef double (*func)(int, const double*, int, const double*, int);
int main()
{
    double a[3];
    double b[3];
    std::cin >> a[0];
    std::cin >> a[1];
    std::cin >> a[2];
    std::cin >> b[0];
    std::cin >> b[1];
    std::cin >> b[2];

    func cblas_ddot;
    HINSTANCE rt = LoadLibrary(TEXT("mkl_rt.dll"));
    cblas_ddot = (func)GetProcAddress(rt, "cblas_ddot");
    double res_rt = cblas_ddot(3, a, 1, b, 1);
    std::cout.precision(25);
    std::cout << res_rt;
}

让我们尝试使用 Visual Studio 的汇编调试器在每台机器上运行它，从 i7（机器 1）/仅支持 AVX2 的机器开始；在这里，在每种情况下，我们都注意到指令修改的所有 YMM 寄存器；例如，YMM4和分别用和YMM5的值进行初始化，在 之后，包含两个数组的元素乘积，在 之后， 的下半部分包含结果：abvfmadd231pdYMM3vaddsdYMM5

vmaskmovpd  ymm4,ymm5,ymmword ptr [rbx]
YMM4 = 0000000000000000-BFEF9C656BB84218-BFB3E64ABC939CC1-3FC13BC946A68994 
vmaskmovpd  ymm5,ymm5,ymmword ptr [r9]
YMM5 = 0000000000000000-BFE026E9B3AD5464-BFA3057691D85EDE-BFEB9951B813250D 
vfmadd231pd ymm3,ymm5,ymm4
YMM3 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC
vaddpd      ymm1,ymm3,ymm1  
YMM1 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC 
vaddpd      ymm0,ymm2,ymm0  
vaddpd      ymm2,ymm1,ymm0
YMM2 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC 
vhaddpd     ymm3,ymm2,ymm2
YMM3 = 3FDFE946951928C9-3FDFE946951928C9-BFBCFCC53F6313B2-BFBCFCC53F6313B2 
vperm2f128  ymm4,ymm3,ymm3,1
YMM4 = BFBCFCC53F6313B2-BFBCFCC53F6313B2-3FDFE946951928C9-3FDFE946951928C9 
vaddsd      xmm5,xmm3,xmm4
YMM5 = 0000000000000000-0000000000000000-BFBCFCC53F6313B2-3FD8AA15454063DC
vmovsd      qword ptr [rsp+90h],xmm5

在支持 AVX-512 的机器 2 上进行相同的实验，结果如下（我们只给出了 ZMM 寄存器的下半部分）：

vmovupd     zmm5{k1}{z},zmmword ptr [r12]
ZMM5 = 0000000000000000-BFEF9C656BB84218-BFB3E64ABC939CC1-3FC13BC946A68994 
vmovupd     zmm4{k1}{z},zmmword ptr [r9]  
ZMM4 = 0000000000000000-BFE026E9B3AD5464-BFA3057691D85EDE-BFEB9951B813250D 
vfmadd231pd zmm3,zmm4,zmm5
ZMM3 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC
vaddpd      zmm17,zmm1,zmm0  
mov         eax,0F0h  
kmovw       k1,eax  
vaddpd      zmm16,zmm3,zmm2
ZMM16= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC
vaddpd      zmm19,zmm16,zmm17 
ZMM19= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC 
mov         eax,0Ch  
kmovw       k2,eax  
vcompresspd zmm18{k1}{z},zmm19  
vaddpd      zmm21,zmm18,zmm19 
ZMM21= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFEC
vcompresspd zmm20{k2}{z},zmm21
ZMM20= 0000000000000000-0000000000000000-0000000000000000-3FDFE946951928C9
vaddpd      zmm0,zmm20,zmm21
ZMM0 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-3FD87AC4BCE238CE 
vhaddpd     xmm1,xmm0,xmm0
ZMM1 = 0000000000000000-0000000000000000-3FD8AA15454063DD-3FD8AA15454063DD 
vmovsd      qword ptr [rsp+88h],xmm1

比较两者，我们首先注意到差异是一个位，3FD8AA15454063DCvs 3FD8AA15454063DD.，但我们现在也看到它是如何出现的：在 AVX2 的情况下，我们对对应于向量的第 0 和第 1 个条目执行水平相加，而在 AVX-512 案例中，我们使用的是第 0 个和第 2 个条目。也就是说，这种差异似乎只是归结为您通过天真计算得到的结果v[0]*u[0] + v[2]*u[2] + v[1]*u[1]与v[0]*u[0] + v[1]*u[1] + v[2]*u[2]. 事实上，比较两者，我们发现完全相同的差异：

In [34]: '%.25f' % (v[0]*u[0] + v[2]*u[2] + v[1]*u[1])
Out[34]: '0.3853810478481685675156143'

In [35]: '%.25f' % (v[0]*u[0] + v[1]*u[1] + v[2]*u[2])
Out[35]: '0.3853810478481685120044631'

推荐阅读