scipy - 消除两台机器之间ddot的差异

我目前有两台机器，它们在两个向量上为np.dot的实例产生不同的输出。无需深入研究从NumPy到BLAS的抽象层次，我就可以再现scipy.linalg.blas.ddot中的差异，因此我假设对BLAS中差异的解释也解释了NumPy中的差异。具体而言，请考虑以下示例：import numpy as npfrom scipy.linalg.blas import ddotu = 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 2a 0.3853810478481685120044631 0.3853810478481685675156143b 0.3853810478481685120044631 0.3853810478481685120044631c 0.3853810478481685120044631 0.3853810478481685675156143同样，以下Cython也会引起相同的差异：cimport scipy.linalg.cython_blascimport numpy as npimport numpy as npcdef 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）。在这两种情况下，我都只设置了numpy，scipy，cython及其依赖项而没有设置新的conda环境。我已经在环境上运行了校验和，以确保最终包含的二进制文件相符并验证了np.__config__.show()的输出是否匹配。同样，我检查了mkl.get_version_string()的输出是否在两台计算机上一致。这使我认为问题可能出在硬件差异上。我没有研究最终执行了哪些指令（缺少在Windows / MSVC上调试Cython代码的直接方法），但是我checked这两台机器都支持AVX2 / FMA，这似乎可能是AVX2 / FMA的一个来源。差异。另一方面，我确实发现这两台机器支持不同的指令集。具体来说 Machine 1 (i7) Machine 2 (Xeon)AVX Y YAVX2 Y YAVX512CD N YAVX512ER N NAVX512F N YAVX512PF N NFMA Y Y但是，我不知道一种确定它本身是否足以解释差异的好方法，或者它是否是红色鲱鱼（？）所以我的问题变成： 从上面开始，有什么自然的步骤可以尝试找出导致差异的原因？是组装时间，还是还有其他更明显的东西？ 最佳答案 给出了对该问题的出色评论，显然，受支持的指令集之间的差异最终是罪魁祸首，的确，我们可以在运行Cython脚本时使用ListDLLs来发现MKL基于这两种情况加载了不同的库。对于i7（机器1）：>listdlls64 python.exe | wsl grep mkl0x00000000b9ff0000 0xe7e000 [...]\miniconda3\envs\npfloattest\Library\bin\mkl_rt.dll0x00000000b80e0000 0x1f05000 [...]\miniconda3\envs\npfloattest\Library\bin\mkl_intel_thread.dll0x00000000b3b40000 0x43ba000 [...]\miniconda3\envs\npfloattest\Library\bin\mkl_core.dll0x00000000b0e50000 0x2ce5000 [...]\miniconda3\envs\npfloattest\Library\bin\mkl_avx2.dll0x00000000b01f0000 0xc58000 [...]\miniconda3\envs\npfloattest\Library\bin\mkl_vml_avx2.dll0x00000000f88c0000 0x7000 [...]\miniconda3\envs\npfloattest\lib\site-packages\mkl\_mklinit.cp37-win_amd64.pyd0x00000000afce0000 0x22000 [...]\miniconda3\envs\npfloattest\lib\site-packages\mkl\_py_mkl_service.cp37-win_amd64.pyd对于Xeon（机器2）：0x0000000057ec0000 0xe7e000 [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_rt.dll0x0000000055fb0000 0x1f05000 [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_intel_thread.dll0x0000000051bf0000 0x43ba000 [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_core.dll0x000000004e1a0000 0x3a4a000 [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_avx512.dll0x000000005c6c0000 0xc03000 [...]\Miniconda3\envs\npfloattest\Library\bin\mkl_vml_avx512.dll0x0000000079a70000 0x7000 [...]\Miniconda3\envs\npfloattest\lib\site-packages\mkl\_mklinit.cp37-win_amd64.pyd0x000000005e830000 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分别用a和b的值初始化，在vfmadd231pd之后，YMM3包含两个数组的元素乘积，而，vaddsd的下部包含结果：vmaskmovpd ymm4,ymm5,ymmword ptr [rbx]YMM4 = 0000000000000000-BFEF9C656BB84218-BFB3E64ABC939CC1-3FC13BC946A68994vmaskmovpd ymm5,ymm5,ymmword ptr [r9]YMM5 = 0000000000000000-BFE026E9B3AD5464-BFA3057691D85EDE-BFEB9951B813250Dvfmadd231pd ymm3,ymm5,ymm4YMM3 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvaddpd ymm1,ymm3,ymm1YMM1 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvaddpd ymm0,ymm2,ymm0vaddpd ymm2,ymm1,ymm0YMM2 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvhaddpd ymm3,ymm2,ymm2YMM3 = 3FDFE946951928C9-3FDFE946951928C9-BFBCFCC53F6313B2-BFBCFCC53F6313B2vperm2f128 ymm4,ymm3,ymm3,1YMM4 = BFBCFCC53F6313B2-BFBCFCC53F6313B2-3FDFE946951928C9-3FDFE946951928C9vaddsd xmm5,xmm3,xmm4YMM5 = 0000000000000000-0000000000000000-BFBCFCC53F6313B2-3FD8AA15454063DCvmovsd qword ptr [rsp+90h],xmm5在机器2上进行的同一实验（一个支持AVX-512的实验）给出了以下结果（此处仅给出了ZMM寄存器的下半部分）：vmovupd zmm5{k1}{z},zmmword ptr [r12]ZMM5 = 0000000000000000-BFEF9C656BB84218-BFB3E64ABC939CC1-3FC13BC946A68994vmovupd zmm4{k1}{z},zmmword ptr [r9]ZMM4 = 0000000000000000-BFE026E9B3AD5464-BFA3057691D85EDE-BFEB9951B813250Dvfmadd231pd zmm3,zmm4,zmm5ZMM3 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvaddpd zmm17,zmm1,zmm0mov eax,0F0hkmovw k1,eaxvaddpd zmm16,zmm3,zmm2ZMM16= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvaddpd zmm19,zmm16,zmm17ZMM19= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECmov eax,0Chkmovw k2,eaxvcompresspd zmm18{k1}{z},zmm19vaddpd zmm21,zmm18,zmm19ZMM21= 0000000000000000-3FDFE946951928C9-3F67A8442F158742-BFBDBA0760DBBFECvcompresspd zmm20{k2}{z},zmm21ZMM20= 0000000000000000-0000000000000000-0000000000000000-3FDFE946951928C9vaddpd zmm0,zmm20,zmm21ZMM0 = 0000000000000000-3FDFE946951928C9-3F67A8442F158742-3FD87AC4BCE238CEvhaddpd xmm1,xmm0,xmm0ZMM1 = 0000000000000000-0000000000000000-3FD8AA15454063DD-3FD8AA15454063DDvmovsd qword ptr [rsp+88h],xmm1比较两者，我们首先注意到差异只是一个位，YMM5 vs. 3FD8AA15454063DC，但是现在我们也看到了它的产生方式：在AVX2情况下，我们对与0th和0th相对应的值执行水平加法。向量的第1个条目，而在AVX-512情况下，我们使用第0个和第2个条目。就是说，差异似乎可以归结为简单地计算3FD8AA15454063DD和v[0]*u[0] + v[2]*u[2] + v[1]*u[1]所得到的差异。确实，比较两者，我们发现完全相同的差异：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'