def block_outer_sum0(x1, x2):
    """
    Multiply x1 (J1 x M) and x2 (J2xM) and extend the dimensions to 3D (J1xJ2xM)
    """
    J1, M = x1.shape
    J2, M = x2.shape
    xx1 = x1.reshape((J1, 1, M), order='C')
    xx2 = x2.reshape((1, J2, M), order='C')
    y = xx1 + xx2
    return y
