Python code inconsistency on complex multiplication in MPI4py

Assume a Python MPI program where a master node sends a pair of complex matrices to each worker node and the worker node is supposed to compute their product (conventional matrix product). The input matrices are constructed at the master node according to some algorithm which there is no need to explain. Now imagine for simplicity that we have only 2 MPI processes, one master and one worker. I have created two versions of this program for this case. The first one constructs two complex numbers (1-by-1 matrices for simplicity) and sends them to the worker to compute the product. This program is like a skeleton for what I am trying to do with multiple workers. In the second program, I have omitted the algorithm and have just hard-coded these two complex numbers into the code. The programs are supposed to give the same product shown here:

a = 28534314.10478439+28534314.10478436j

b = -1.39818115e+09+1.39818115e+09j

a*b = -7.97922802e+16+48j

This has been checked in Matlab. Instead, the first program does not work and the worker gives a*b = -7.97922801e+16+28534416.j while the second program works correctly. Please note that the data is transmitted correctly from the master to the worker and the data structures are the same in both cases (see the print() functions).

The first (wrong) program is:

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(np.matrix(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(np.matrix(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)

    #Test
    # print(type(a[0]))
    # print(a[0].dtype)
    print("Sender sent:  ")
    print(a[0])
    print(b[0])

else:

    print("I am receiver")

    A = np.zeros_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
    B = np.zeros_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)

The second (correct) program is:

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').astype(np.complex128)
    b = np.matrix('-1.39818115e+09+1.39818115e+09j').astype(np.complex128)

    #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)

    #Test
    # print(type(a[0]))
    # print(a[0].dtype)
    print("Sender sent:  ")
    print(a[0])
    print(b[0])

else:

    print("I am receiver")

    A = np.zeros_like(np.matrix([[0]*(1) for i in range(1)])).astype(np.complex128)
    B = np.zeros_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)

I am using MPI4py 3.0.0. along with Python 2.7.14 and the kernel of Open MPI 2.1.2. I have been straggling with this problem for a whole day and still cannot figure out what's going on. I have tried numerous initializations like np.zeros(), np.zeros_like(), np.empty_like() as well as both np.array and np.matrix and functions np.dot(), np.matmul() and the operator *. I have also tried doing the multiplication using sympy or Python objects.

Finally, I think that the problem is always with the imaginary part of the product based on other examples I tried.

This seems to be a major issue if I am not missing something.

Comments (3)