def split_matrix(rdm):
    N = rdm.rdm.shape[0]
    R = int(np.sqrt(N))
    nn = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    ne = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    ee = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    for i in range(N):
        (p, r) = ((i // R), (i % R))
        for j in range(N):
            (q, s) = ((j // R), (j % R))
            ind = tuple([p, q, r, s])
            if (len(set(ind)) == 2):
                nn[(i, j)] = rdm.rdm[(i, j)]
            elif (len(set(ind)) == 3):
                ne[(i, j)] = rdm.rdm[(i, j)]
            elif (len(set(ind)) == 4):
                ee[(i, j)] = rdm.rdm[(i, j)]
    return (nn, ne, ee)
def split_matrix(rdm):
    N = rdm.rdm.shape[0]
    R = int(np.sqrt(N))
    nn = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    ne = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    ee = np.zeros(rdm.rdm.shape, dtype=np.complex_)
    for i in range(N):
        (p, r) = ((i // R), (i % R))
        for j in range(N):
            (q, s) = ((j // R), (j % R))
            ind = tuple([p, q, r, s])
            if (len(set(ind)) == 2):
                nn[(i, j)] = rdm.rdm[(i, j)]
            elif (len(set(ind)) == 3):
                ne[(i, j)] = rdm.rdm[(i, j)]
            elif (len(set(ind)) == 4):
                ee[(i, j)] = rdm.rdm[(i, j)]
    return (nn, ne, ee)