def peo(A, p):
    """
    Checks whether an ordering is a perfect elmimination order.

    Returns `True` if the permutation :math:`p` is a perfect elimination order
    for a Cholesky factorization :math:`PAP^T = LL^T`. Only the lower
    triangular part of :math:`A` is accessed.
    
    :param A:   :py:class:`spmatrix`
    
    :param p:   :py:class:`matrix` or :class:`list` of length `A.size[0]`
    """

    n = A.size[0]
    assert type(A) == spmatrix, "A must be a sparse matrix"
    assert A.size[1] == n, "A must be a square matrix"
    assert len(p) == n, "length of p must be equal to the order of A" 
    if isinstance(p, list): p = matrix(p)
    
    As = symmetrize(A)
    cp,ri,_ = As.CCS

    # compute inverse permutation array
    ip = matrix(0,(n,1))
    ip[p] = matrix(range(n),(n,1))

    # test set inclusion
    for k in range(n):
        v = p[k]  # next vertex to be eliminated

        # indices of neighbors that correspond to strictly lower triangular elements in reordered pattern
        r = set([rj for rj in ri[cp[v]:cp[v+1]] if ip[rj] > k])  

        for rj in r:
            if not r.issubset(set(ri[cp[rj]:cp[rj+1]])): return False
            
    return True
Checks whether an ordering is a perfect elmimination order.

    Returns `True` if the permutation :math:`p` is a perfect elimination order
    for a Cholesky factorization :math:`PAP^T = LL^T`. Only the lower
    triangular part of :math:`A` is accessed.
    
    :param A:   :py:class:`spmatrix`
    
    :param p:   :py:class:`matrix` or :class:`list` of length `A.size[0]`