确定两台计算机之间的 ddot 差异

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

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

但是，我不知道有什么好方法来确定这本身是否足以解释差异，或者它是否是红鲱鱼(？

所以我的问题变成了：

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

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

对于 i7(机器 1(：

>listdlls64 python.exe | wsl grep mkl
0x00000000b9ff0000  0xe7e000   [...]miniconda3envsnpfloattestLibrarybinmkl_rt.dll
0x00000000b80e0000  0x1f05000  [...]miniconda3envsnpfloattestLibrarybinmkl_intel_thread.dll
0x00000000b3b40000  0x43ba000  [...]miniconda3envsnpfloattestLibrarybinmkl_core.dll
0x00000000b0e50000  0x2ce5000  [...]miniconda3envsnpfloattestLibrarybinmkl_avx2.dll
0x00000000b01f0000  0xc58000   [...]miniconda3envsnpfloattestLibrarybinmkl_vml_avx2.dll
0x00000000f88c0000  0x7000     [...]miniconda3envsnpfloattestlibsite-packagesmkl_mklinit.cp37-win_amd64.pyd
0x00000000afce0000  0x22000    [...]miniconda3envsnpfloattestlibsite-packagesmkl_py_mkl_service.cp37-win_amd64.pyd

对于至强(机器 2(：

0x0000000057ec0000  0xe7e000   [...]Miniconda3envsnpfloattestLibrarybinmkl_rt.dll
0x0000000055fb0000  0x1f05000  [...]Miniconda3envsnpfloattestLibrarybinmkl_intel_thread.dll
0x0000000051bf0000  0x43ba000  [...]Miniconda3envsnpfloattestLibrarybinmkl_core.dll
0x000000004e1a0000  0x3a4a000  [...]Miniconda3envsnpfloattestLibrarybinmkl_avx512.dll
0x000000005c6c0000  0xc03000   [...]Miniconda3envsnpfloattestLibrarybinmkl_vml_avx512.dll
0x0000000079a70000  0x7000     [...]Miniconda3envsnpfloattestlibsite-packagesmkl_mklinit.cp37-win_amd64.pyd
0x000000005e830000  0x22000    [...]Miniconda3envsnpfloattestlibsite-packagesmkl_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之后，YMM5的下半部分包含结果：

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

在机器 2 上(支持 AVX-512 的机器(上的相同实验给出了以下结果(我们只给出了 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

比较两者，我们首先注意到差异是一位，3FD8AA15454063DC与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'

如果您需要逐位相等的答案：您是否尝试过从@bg2b(https://software.intel.com/en-us/mkl-linux-developer-guide-instruction-set-specific-dispatching-on-intel-architectures(发布的链接中MKL_ENABLE_INSTRUCTIONS变量？如果您导入 MKL 库然后调用mkl.enable_instructions则可能为时已晚。
在双精度(DP(世界中：相对差值为-1.4404224470807435333001684021155e-16(绝对差值为-5.55111512e-17(，小于C++和Python DP机epsilon(https://en.wikipedia.org/wiki/Machine_epsilon(。因此，结果与Python的正确性相等。

干杯弗拉基米尔

相关内容

最新更新

热门标签：