Python MPI 中复数乘法的代码不一致
Python code inconsistency on complex multiplication in MPI
假设一个 Python MPI 程序,其中主节点向每个工作节点发送一对复杂矩阵,并且工作节点应该计算它们的乘积(常规矩阵乘积)。输入矩阵是根据某种无需解释的算法在主节点处构建的。现在为简单起见,假设我们只有 2 个 MPI 进程,一个主进程和一个工作进程。我为此案例创建了该程序的两个版本。第一个构造两个复数(为简单起见,为 1×1 矩阵)并将它们发送给 worker 以计算乘积。这个程序就像是我试图与多个工人一起做的事情的骨架。在第二个程序中,我省略了算法,只是将这两个复数硬编码到代码中。这些程序应该提供与此处所示相同的产品:
a = 28534314.10478439+28534314.10478436j
b = -1.39818115e+09+1.39818115e+09j
a*b = -7.97922802e+16+48j
这已经在 Matlab 中检查过了。相反,第一个程序不工作,工作人员给出 a*b = -7.97922801e+16+28534416.j
而第二个程序正常工作。请注意,在这两种情况下,数据都从 master 正确传输到 worker(参见 print()
函数)。第一个(错误的)程序是:
from mpi4py import MPI
import numpy as np
N = 1
ell = 9
s_cod = 7
var = np.array([np.exp(1j*2*np.pi*1/8)])
comm = MPI.COMM_WORLD
if comm.rank == 0:
print("I am sender")
#Construction algorithm, explanation skipped
A=np.matrix('1 0; 0 1')
B=np.matrix('1 0; 0 1')
Ah=np.split(A,2)
Bh=np.split(B,2)
Ahv = []
Bhv = []
for i in range(2):
Ahv.append(np.split(Ah[i], 2, axis=1))
Bhv.append(np.split(Bh[i], 2, axis=1))
a = []
b = []
for i in range(N):
a.append(Ahv[0][0]*(pow(s_cod*var[i], ell)) + Ahv[1][0] + Ahv[0][1]*(pow(s_cod*var[i], ell+1)) + Ahv[1][1]*s_cod*var[i])
b.append(Bhv[0][0] + Bhv[1][0]*(pow(s_cod*var[i], ell)) + Bhv[0][1]*(pow(s_cod*var[i], 2)) + Bhv[1][1]*(pow(s_cod*var[i], ell+2)))
#Send message with a predefined tag, like 15 and 16, to each receiver
for i in range(N):
comm.Isend([a[i],MPI.COMPLEX], dest=i+1, tag=15)
comm.Isend([b[i],MPI.COMPLEX], dest=i+1, tag=16)
print("Sender sent: ")
print(a[0])
print(b[0])
else:
print("I am receiver")
A = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
B = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
#Receive message with tags 15, 16 from rank 0
rA = comm.Irecv(A, source=0, tag=15)
rB = comm.Irecv(B, source=0, tag=16)
rA.wait()
rB.wait()
C = np.dot(A, B)
print("Receiver received: ")
print(A)
print(B)
print("Receiver computed: ")
print(C)
第二个(正确的)程序是:
from mpi4py import MPI
import numpy as np
comm = MPI.COMM_WORLD
if comm.rank == 0:
print("I am sender")
a = np.matrix('28534314.10478439+28534314.10478436j')
b = np.matrix('-1.39818115e+09+1.39818115e+09j')
#Send message with a predefined tag, like 15 and 16, to rank 1
comm.Isend([a, MPI.COMPLEX], dest=1, tag=15)
comm.Isend([b, MPI.COMPLEX], dest=1, tag=16)
print("Sender sent: ")
print(a[0])
print(b[0])
else:
print("I am receiver")
A = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
B = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
#Receive message with tags 15, 16 from rank 0
rA = comm.Irecv(A, source=0, tag=15)
rB = comm.Irecv(B, source=0, tag=16)
rA.wait()
rB.wait()
C = np.dot(A, B)
print("Receiver received: ")
print(A)
print(B)
print("Receiver computed: ")
print(C)
我正在使用 MPI4py 3.0.0。以及 Python 2.7.14 和 Open MPI 2.1.2 的内核。我一整天都在为这个问题苦苦挣扎,但仍然无法弄清楚发生了什么。我尝试了很多初始化,例如 np.zeros()
、np.zeros_like()
、np.empty_like()
以及 np.array
和 np.matrix
以及函数 np.dot()
、np.matmul()
和运算符 *
。最后,我认为问题总是出在基于我尝试过的其他示例的产品的虚部。有什么建议吗?
这与 MPI 完全无关。
np.set_printoptions(precision=15)
确认计算出的 a
和 b
实际上与您输入 "correct" 版本的不同。
我不确定结果的基本事实是什么。计算过程中可能会出现影响增大的舍入误差。在点积过程中差异显着,因为在 "correct" 版本中 b
的 real/imaginary 部分的绝对值是相等的,在计算版本中它们只是相对接近但是有一个显着绝对差。
假设一个 Python MPI 程序,其中主节点向每个工作节点发送一对复杂矩阵,并且工作节点应该计算它们的乘积(常规矩阵乘积)。输入矩阵是根据某种无需解释的算法在主节点处构建的。现在为简单起见,假设我们只有 2 个 MPI 进程,一个主进程和一个工作进程。我为此案例创建了该程序的两个版本。第一个构造两个复数(为简单起见,为 1×1 矩阵)并将它们发送给 worker 以计算乘积。这个程序就像是我试图与多个工人一起做的事情的骨架。在第二个程序中,我省略了算法,只是将这两个复数硬编码到代码中。这些程序应该提供与此处所示相同的产品:
a = 28534314.10478439+28534314.10478436j
b = -1.39818115e+09+1.39818115e+09j
a*b = -7.97922802e+16+48j
这已经在 Matlab 中检查过了。相反,第一个程序不工作,工作人员给出 a*b = -7.97922801e+16+28534416.j
而第二个程序正常工作。请注意,在这两种情况下,数据都从 master 正确传输到 worker(参见 print()
函数)。第一个(错误的)程序是:
from mpi4py import MPI
import numpy as np
N = 1
ell = 9
s_cod = 7
var = np.array([np.exp(1j*2*np.pi*1/8)])
comm = MPI.COMM_WORLD
if comm.rank == 0:
print("I am sender")
#Construction algorithm, explanation skipped
A=np.matrix('1 0; 0 1')
B=np.matrix('1 0; 0 1')
Ah=np.split(A,2)
Bh=np.split(B,2)
Ahv = []
Bhv = []
for i in range(2):
Ahv.append(np.split(Ah[i], 2, axis=1))
Bhv.append(np.split(Bh[i], 2, axis=1))
a = []
b = []
for i in range(N):
a.append(Ahv[0][0]*(pow(s_cod*var[i], ell)) + Ahv[1][0] + Ahv[0][1]*(pow(s_cod*var[i], ell+1)) + Ahv[1][1]*s_cod*var[i])
b.append(Bhv[0][0] + Bhv[1][0]*(pow(s_cod*var[i], ell)) + Bhv[0][1]*(pow(s_cod*var[i], 2)) + Bhv[1][1]*(pow(s_cod*var[i], ell+2)))
#Send message with a predefined tag, like 15 and 16, to each receiver
for i in range(N):
comm.Isend([a[i],MPI.COMPLEX], dest=i+1, tag=15)
comm.Isend([b[i],MPI.COMPLEX], dest=i+1, tag=16)
print("Sender sent: ")
print(a[0])
print(b[0])
else:
print("I am receiver")
A = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
B = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
#Receive message with tags 15, 16 from rank 0
rA = comm.Irecv(A, source=0, tag=15)
rB = comm.Irecv(B, source=0, tag=16)
rA.wait()
rB.wait()
C = np.dot(A, B)
print("Receiver received: ")
print(A)
print(B)
print("Receiver computed: ")
print(C)
第二个(正确的)程序是:
from mpi4py import MPI
import numpy as np
comm = MPI.COMM_WORLD
if comm.rank == 0:
print("I am sender")
a = np.matrix('28534314.10478439+28534314.10478436j')
b = np.matrix('-1.39818115e+09+1.39818115e+09j')
#Send message with a predefined tag, like 15 and 16, to rank 1
comm.Isend([a, MPI.COMPLEX], dest=1, tag=15)
comm.Isend([b, MPI.COMPLEX], dest=1, tag=16)
print("Sender sent: ")
print(a[0])
print(b[0])
else:
print("I am receiver")
A = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
B = np.empty_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
#Receive message with tags 15, 16 from rank 0
rA = comm.Irecv(A, source=0, tag=15)
rB = comm.Irecv(B, source=0, tag=16)
rA.wait()
rB.wait()
C = np.dot(A, B)
print("Receiver received: ")
print(A)
print(B)
print("Receiver computed: ")
print(C)
我正在使用 MPI4py 3.0.0。以及 Python 2.7.14 和 Open MPI 2.1.2 的内核。我一整天都在为这个问题苦苦挣扎,但仍然无法弄清楚发生了什么。我尝试了很多初始化,例如 np.zeros()
、np.zeros_like()
、np.empty_like()
以及 np.array
和 np.matrix
以及函数 np.dot()
、np.matmul()
和运算符 *
。最后,我认为问题总是出在基于我尝试过的其他示例的产品的虚部。有什么建议吗?
这与 MPI 完全无关。
np.set_printoptions(precision=15)
确认计算出的 a
和 b
实际上与您输入 "correct" 版本的不同。
我不确定结果的基本事实是什么。计算过程中可能会出现影响增大的舍入误差。在点积过程中差异显着,因为在 "correct" 版本中 b
的 real/imaginary 部分的绝对值是相等的,在计算版本中它们只是相对接近但是有一个显着绝对差。