Error occurred when calling function dot with transA/B='H' option

[zuoshifan]

When calling function dot with transA/B='H' option, strange error occurs like

--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD 
with errorcode -2.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
[node2:09438] 3 more processes have sent help message help-mpi-api.txt / mpi-abort
[node2:09438] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

Here is my test script

import numpy as np
import scipy.linalg as la

from mpi4py import MPI

import scalapy
from scalapy import core
import scalapy.routines as rt

comm = MPI.COMM_WORLD

rank = comm.rank
size = comm.size

if size != 4:
    raise Exception("Test needs 4 processes.")

core.initmpi([2, 2], block_shape=[3, 3])

allclose = lambda a, b: np.allclose(a, b, rtol=1e-4, atol=1e-6)

def test_dot_D():
    ns = 5

    gA = np.random.standard_normal((ns, ns)).astype(np.float64)
    gA = np.asfortranarray(gA)
    gB = np.random.standard_normal((ns, ns)).astype(np.float64)
    gB = np.asfortranarray(gB)

    dA = core.DistributedMatrix.from_global_array(gA, rank=0)
    dB = core.DistributedMatrix.from_global_array(gB, rank=0)

    dAdB = rt.dot(dA, dB, transA='N', transB='N')
    dAdB = dAdB.to_global_array(rank=0)
    dAdBT = rt.dot(dA, dB, transA='N', transB='T')
    dAdBT = dAdBT.to_global_array(rank=0)
    dATdB = rt.dot(dA, dB, transA='T', transB='N')
    dATdB = dATdB.to_global_array(rank=0)
    dATdBT = rt.dot(dA, dB, transA='T', transB='T')
    dATdBT = dATdBT.to_global_array(rank=0)

    if rank == 0:
        gAgB = np.dot(gA, gB)
        gAgBT = np.dot(gA, gB.T)
        gATgB = np.dot(gA.T, gB)
        gATgBT = np.dot(gA.T, gB.T)

        assert allclose(dAdB, gAgB)
        assert allclose(dAdBT, gAgBT)
        assert allclose(dATdB, gATgB)
        assert allclose(dATdBT, gATgBT)

def test_dot_Z():
    ns = 5

    gA = np.random.standard_normal((ns, ns)).astype(np.float64)
    gA = gA + 1.0J * np.random.standard_normal((ns, ns)).astype(np.float64)
    gA = np.asfortranarray(gA)
    gB = np.random.standard_normal((ns, ns)).astype(np.float64)
    gB = gB + 1.0J * np.random.standard_normal((ns, ns)).astype(np.float64)
    gB = np.asfortranarray(gB)

    dA = core.DistributedMatrix.from_global_array(gA, rank=0)
    dB = core.DistributedMatrix.from_global_array(gB, rank=0)

    dAdB = rt.dot(dA, dB, transA='N', transB='N')
    dAdB = dAdB.to_global_array(rank=0)
    dAdBH = rt.dot(dA, dB, transA='N', transB='H')
    dAdBH = dAdBH.to_global_array(rank=0)
    dAHdB = rt.dot(dA, dB, transA='H', transB='N')
    dAHdB = dAHdB.to_global_array(rank=0)
    dAHdBH = rt.dot(dA, dB, transA='H', transB='H')
    dAHdBH = dAHdBH.to_global_array(rank=0)

    if rank == 0:
        gAgB = np.dot(gA, gB)
        gAgBH = np.dot(gA, gB.T.conj())
        gAHgB = np.dot(gA.T.conj(), gB)
        gAHgBH = np.dot(gA.T.conj(), gB.T.conj())

        assert allclose(dAdB, gAgB)
        assert allclose(dAdBH, gAgBH)
        assert allclose(dAHdB, gAHgB)
        assert allclose(dAHdBH, gAHgBH)

if __name__ == '__main__':
    test_dot_D()
    test_dot_Z()

[zuoshifan]

It seems that the option transA/B='C' is for conjugate transpose (It is the ScaLapack convention, too, see http://www.netlib.org/scalapack/slug/node95.html), so the docstring

transA, transB : ['N', 'T', 'H', 'C']
        Whether we should use a transpose, rather than A or B themselves.
        Either, do nothing ('N'), normal transpose ('T'), Hermitian transpose
        ('H'), or complex conjugation only ('C').

is deceptive.

The test code snippet is

dAdBC = rt.dot(dA, dB, transA='N', transB='C')
dAdBC = dAdBC.to_global_array(rank=0)
...
if rank == 0:
        gAgBC = np.dot(gA, gB.T.conj())

        assert allclose(dAdBC, gAgBC)

[jrs65]

Thanks @zuoshifan for pointing this one out. Do you mind if I take your test case from above, and use it to replace the test suite in test_dot.py.

[zuoshifan]

Sure, no problem.