def _P(m):
    """Returns m^2 x m^2 permutation matrix that swaps rows i and j where
    j = 1 + m((i - 1) mod m) + (i - 1) div m, for i = 1 .. m^2
    """
    P = np.zeros((m**2,m**2), dtype=np.int64)
    for i in range(1, m**2 + 1):
        j = 1 + m*((i - 1) % m) + (i - 1)//m
        P[i-1, j-1] = 1
    return P