Maps a symbolic factorization to a block-diagonal structure with
    coupling constraints.

    :param symb:               :py:class:`symbolic`
    :param coupling:           optional 
    :return dims:              list of block dimensions
    :return sparse_to_block:   dictionary
    :return constraints:       list of coupling constraints

def symb_to_block(symb, coupling = 'full'):
    """
    Maps a symbolic factorization to a block-diagonal structure with
    coupling constraints.

    :param symb:               :py:class:`symbolic`
    :param coupling:           optional 
    :return dims:              list of block dimensions
    :return sparse_to_block:   dictionary
    :return constraints:       list of coupling constraints
    """

    n = len(symb.snode)             # order of block
    Ncliques = len(symb.snpar)      # number of cliques

    # compute clique orders
    dims = [symb.sncolptr[j+1]-symb.sncolptr[j] for j in range(Ncliques)]

    # compute offsets in block-diagonal structure
    offsets = [0]
    for i in range(Ncliques): offsets.append(offsets[-1] + dims[i]**2)

    constraints = []        # list of coupling constraints
    sparse_to_block = {}    # conversion dictionary

    for k in range(Ncliques):
        # map nonzeros in {Jk,Nk} part of clique k to block-diagonal structure
        nodes = symb.snode[symb.snptr[k]:symb.snptr[k+1]]
        rows = symb.snrowidx[symb.sncolptr[k]:symb.sncolptr[k+1]]
        nk = len(nodes)    # number of nodes in supernode
        wk = len(rows)     # number of nodes in clique
        for j in range(nk):
            for i in range(j,wk):
                if i == j:
                    sparse_to_block[nodes[j]*n + rows[i]] = (offsets[k] + j*wk + i,)
                else:
                    sparse_to_block[nodes[j]*n + rows[i]] =(offsets[k] + j*wk + i, offsets[k] + i*wk + j)

        # add coupling constraints to list of constraints
        if symb.snpar[k] == k: continue   # skip if supernode k is a root supernode
        p = symb.snpar[k]
        np = len(symb.snode[symb.snptr[p]:symb.snptr[p+1]])
        wp = symb.sncolptr[p+1] - symb.sncolptr[p]
        ri = symb.relidx[symb.relptr[k]:symb.relptr[k+1]]

        if type(coupling) is spmatrix:           
            tmp = coupling[rows[nk:],rows[nk:]]
            for i,j in zip(tmp.I,tmp.J):
                if j == i:
                    constraints.append((offsets[k] + (j+nk)*wk + i+nk,
                                        offsets[p] + ri[j]*wp + ri[i]))
                else:
                    constraints.append((offsets[k] + (j+nk)*wk + i+nk,
                                        offsets[p] + ri[j]*wp + ri[i],
                                        offsets[k] + (i+nk)*wk + j+nk,
                                        offsets[p] + ri[i]*wp + ri[j]))   
                     
        elif coupling == 'full':
            for j in range(len(ri)):
                for i in range(j,len(ri)):
                    if j == i:
                        constraints.append((offsets[k] + (j+nk)*wk + i+nk,
                                            offsets[p] + ri[j]*wp + ri[i]))
                    else:
                        constraints.append((offsets[k] + (j+nk)*wk + i+nk,
                                            offsets[p] + ri[j]*wp + ri[i],
                                            offsets[k] + (i+nk)*wk + j+nk,
                                            offsets[p] + ri[i]*wp + ri[j]))

    return dims, sparse_to_block, constraints