Find supernodes and supernodal etree.

    ARGUMENTS
    par       parent array

    post      array with post ordering

    colcount  array with column counts

    RETURNS
    snpar     supernodal parent structure 

    flag      integer vector of length n; if flag[i] < 0, then -flag[i]
              is the degree of the supernode with repr. vertex i; if
              flag[i] >= 0, then flag[i] is the repr. vertex to which
              node i belongs.

def pothen_sun(par, post, colcount):
    """
    Find supernodes and supernodal etree.

    ARGUMENTS
    par       parent array

    post      array with post ordering

    colcount  array with column counts

    RETURNS
    snpar     supernodal parent structure 

    flag      integer vector of length n; if flag[i] < 0, then -flag[i]
              is the degree of the supernode with repr. vertex i; if
              flag[i] >= 0, then flag[i] is the repr. vertex to which
              node i belongs.
    """
    
    n = len(par)
    flag = matrix(-1, (n, 1))
    snpar = matrix(-1, (n, 1))
    snodes = n
    ch = {}
    
    for j in post:

        if par[j] in ch: ch[par[j]].append(j)
        else: ch[par[j]] = [j]

        mdeg = colcount[j] - 1

        if par[j] != j:
            if mdeg == colcount[par[j]] and flag[par[j]] == -1:
                # par[j] not assigned to supernode
                snodes -= 1
                if flag[j] < 0:   # j is a repr. vertex
                    flag[par[j]] = j
                    flag[j] -= 1
                else:             # j is not a repr. vertex
                    flag[par[j]] = flag[j]
                    flag[flag[j]] -= 1
        else:
            if flag[j] < 0: snpar[j] = j
            else: snpar[flag[j]] = flag[j]

        if flag[j] < 0: k = j
        else: k = flag[j]

        if j in ch:
            for i in ch[j]:
                if flag[i] < 0: l = i
                else: l = flag[i]
                if not l == k: snpar[l] = k


    repr = matrix([i for i in range(n) if flag[i] < 0])
    deg = matrix([-flag[i] for i in range(n) if flag[i] < 0])

    # renumber etree with number of supernodes
    sn = matrix(-1, (n+1, 1))
    for k, r in enumerate(repr): sn[r] = k
    snpar = sn[snpar[repr]]

    return snpar, flag