The R Project SVN R-packages

/* ========================================================================== */
/* === Partition/cholmod_nesdis ============================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Partition Module.
 * Copyright (C) 2005-2006, Univ. of Florida.  Author: Timothy A. Davis
 * The CHOLMOD/Partition Module is licensed under Version 2.1 of the GNU
 * Lesser General Public License.  See lesser.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* CHOLMOD nested dissection and graph partitioning.
 *
 * cholmod_bisect:
 *
 *  Finds a set of nodes that partitions the graph into two parts.
 *  Compresses the graph first.  Requires METIS.
 *
 * cholmod_nested_dissection:
 *
 *  Nested dissection, using its own compression and connected-commponents
 *  algorithms, an external graph partitioner (METIS), and a constrained
 *  minimum degree ordering algorithm (CCOLAMD or CSYMAMD).  Typically
 *  gives better orderings than METIS_NodeND (about 5% to 10% fewer
 *  nonzeros in L).
 *
 * cholmod_collapse_septree:
 *
 *  Prune the separator tree returned by cholmod_nested_dissection.
 *
 * This file contains several routines private to this file:
 *
 *  partition   compress and partition a graph
 *  clear_flag  clear Common->Flag, but do not modify negative entries
 *  find_components find the connected components of a graph
 *
 * Supports any xtype (pattern, real, complex, or zomplex).
 */

#ifndef NPARTITION

#include "cholmod_internal.h"
#include "cholmod_partition.h"
#include "cholmod_cholesky.h"

/* ========================================================================== */
/* === partition ============================================================ */
/* ========================================================================== */

/* Find a set of nodes that partition a graph.  The graph must be symmetric
 * with no diagonal entries.  To compress the graph first, compress is TRUE
 * and on input Hash [j] holds the hash key for node j, which must be in the
 * range 0 to csize-1. The input graph (Cp, Ci) is destroyed.  Cew is all 1's
 * on input and output.  Cnw [j] > 0 is the initial weight of node j.  On
 * output, Cnw [i] = 0 if node i is absorbed into j and the original weight
 * Cnw [i] is added to Cnw [j].  If compress is FALSE, the graph is not
 * compressed and Cnw and Hash are unmodified.  The partition itself is held in
 * the output array Part of size n.  Part [j] is 0, 1, or 2, depending on
 * whether node j is in the left part of the graph, the right part, or the
 * separator, respectively.  Note that the input graph need not be connected,
 * and the output subgraphs (the three parts) may also be unconnected.
 *
 * Returns the size of the separator, in terms of the sum of the weights of
 * the nodes.  It is guaranteed to be between 1 and the total weight of all
 * the nodes.  If it is of size less than the total weight, then both the left
 * and right parts are guaranteed to be non-empty (this guarantee depends on
 * cholmod_metis_bisector).
 */

static UF_long partition    /* size of separator or -1 if failure */
(
    /* inputs, not modified on output */
#ifndef NDEBUG
    Int csize,      /* upper bound on # of edges in the graph;
             * csize >= MAX (n, nnz(C)) must hold. */
#endif
    int compress,   /* if TRUE the compress the graph first */

    /* input/output */
    Int Hash [ ],   /* Hash [i] = hash >= 0 is the hash function for node
             * i on input.  On output, Hash [i] = FLIP (j) if node
             * i is absorbed into j.  Hash [i] >= 0 if i has not
             * been absorbed. */

    /* input graph, compressed graph of cn nodes on output */
    cholmod_sparse *C,

    /* input/output */
    Int Cnw [ ],    /* size n.  Cnw [j] > 0 is the weight of node j on
             * input.  On output, if node i is absorbed into
             * node j, then Cnw [i] = 0 and the original weight of
             * node i is added to Cnw [j].  The sum of Cnw [0..n-1]
             * is not modified. */

    /* workspace */
    Int Cew [ ],    /* size csize, all 1's on input and output */

    /* more workspace, undefined on input and output */
    Int Cmap [ ],   /* size n (i/i/l) */

    /* output */
    Int Part [ ],   /* size n, Part [j] = 0, 1, or 2. */

    cholmod_common *Common
)
{
    Int n, hash, head, i, j, k, p, pend, ilen, ilast, pi, piend,
    jlen, ok, cn, csep, pdest, nodes_pruned, nz, total_weight, jscattered ;
    Int *Cp, *Ci, *Next, *Hhead ;

#ifndef NDEBUG
    Int cnt, pruned ;
    double work = 0, goodwork = 0 ;
#endif

    /* ---------------------------------------------------------------------- */
    /* quick return for small or empty graphs */
    /* ---------------------------------------------------------------------- */

    n = C->nrow ;
    Cp = C->p ;
    Ci = C->i ;
    nz = Cp [n] ;

    PRINT2 (("Partition start, n "ID" nz "ID"\n", n, nz)) ;

    total_weight = 0 ;
    for (j = 0 ; j < n ; j++)
    {
    ASSERT (Cnw [j] > 0) ;
    total_weight += Cnw [j] ;
    }

    if (n <= 2)
    {
    /* very small graph */
    for (j = 0 ; j < n ; j++)
    {
        Part [j] = 2 ;
    }
    return (total_weight) ;
    }
    else if (nz <= 0)
    {
    /* no edges, this is easy */
    PRINT2 (("diagonal matrix\n")) ;
    k = n/2 ;
    for (j = 0 ; j < k ; j++)
    {
        Part [j] = 0 ;
    }
    for ( ; j < n ; j++)
    {
        Part [j] = 1 ;
    }
    /* ensure the separator is not empty (required by nested dissection) */
    Part [n-1] = 2 ;
    return (Cnw [n-1]) ;
    }

#ifndef NDEBUG
    ASSERT (n > 1 && nz > 0) ;
    PRINT2 (("original graph:\n")) ;
    for (j = 0 ; j < n ; j++)
    {
    PRINT2 ((""ID": ", j)) ;
    for (p = Cp [j] ; p < Cp [j+1] ; p++)
    {
        i = Ci [p] ;
        PRINT3 ((""ID" ", i)) ;
        ASSERT (i >= 0 && i < n && i != j) ;
    }
    PRINT2 (("hash: "ID"\n", Hash [j])) ;
    }
    DEBUG (for (p = 0 ; p < csize ; p++) ASSERT (Cew [p] == 1)) ;
#endif

    nodes_pruned = 0 ;

    if (compress)
    {

    /* ------------------------------------------------------------------ */
    /* get workspace */
    /* ------------------------------------------------------------------ */

    Next = Part ;   /* use Part as workspace for Next [ */
    Hhead = Cew ;   /* use Cew as workspace for Hhead [ */

    /* ------------------------------------------------------------------ */
    /* create the hash buckets */
    /* ------------------------------------------------------------------ */

    for (j = 0 ; j < n ; j++)
    {
        /* get the hash key for node j */
        hash = Hash [j] ;
        ASSERT (hash >= 0 && hash < csize) ;
        head = Hhead [hash] ;
        if (head > EMPTY)
        {
        /* hash bucket for this hash key is empty. */
        head = EMPTY ;
        }
        else
        {
        /* hash bucket for this hash key is not empty.  get old head */
        head = FLIP (head) ;
        ASSERT (head >= 0 && head < n) ;
        }
        /* node j becomes the new head of the hash bucket.  FLIP it so that
         * we can tell the difference between an empty or non-empty hash
         * bucket. */
        Hhead [hash] = FLIP (j) ;
        Next [j] = head ;
        ASSERT (head >= EMPTY && head < n) ;
    }

#ifndef NDEBUG
    for (cnt = 0, k = 0 ; k < n ; k++)
    {
        ASSERT (Hash [k] >= 0 && Hash [k] < csize) ;    /* k is alive */
        hash = Hash [k] ;
        ASSERT (hash >= 0 && hash < csize) ;
        head = Hhead [hash] ;
        ASSERT (head < EMPTY) ; /* hash bucket not empty */
        j = FLIP (head) ;
        ASSERT (j >= 0 && j < n) ;
        if (j == k)
        {
        PRINT2 (("hash "ID": ", hash)) ;
        for ( ; j != EMPTY ; j = Next [j])
        {
            PRINT3 ((" "ID"", j)) ;
            ASSERT (j >= 0 && j < n) ;
            ASSERT (Hash [j] == hash) ;
            cnt++ ;
            ASSERT (cnt <= n) ;
        }
        PRINT2 (("\n")) ;
        }
    }
    ASSERT (cnt == n) ;
#endif

    /* ------------------------------------------------------------------ */
    /* scan the non-empty hash buckets for indistinguishable nodes */
    /* ------------------------------------------------------------------ */

    /* If there are no hash collisions and no compression occurs, this takes
     * O(n) time.  If no hash collisions, but some nodes are removed, this
     * takes time O(n+e) where e is the sum of the degress of the nodes
     * that are removed.  Even with many hash collisions (a rare case),
     * this algorithm has never been observed to perform more than nnz(A)
     * useless work.
     *
     * Cmap is used as workspace to mark nodes of the graph, [
     * for comparing the nonzero patterns of two nodes i and j.
     */

#define Cmap_MARK(i)   Cmap [i] = j
#define Cmap_MARKED(i) (Cmap [i] == j)

    for (i = 0 ; i < n ; i++)
    {
        Cmap [i] = EMPTY ;
    }

    for (k = 0 ; k < n ; k++)
    {
        hash = Hash [k] ;
        ASSERT (hash >= FLIP (n-1) && hash < csize) ;
        if (hash < 0)
        {
        /* node k has already been absorbed into some other node */
        ASSERT (FLIP (Hash [k]) >= 0 && FLIP (Hash [k] < n)) ;
        continue ;
        }
        head = Hhead [hash] ;
        ASSERT (head < EMPTY || head == 1) ;
        if (head == 1)
        {
        /* hash bucket is already empty */
        continue ;
        }
        PRINT2 (("\n--------------------hash "ID":\n", hash)) ;
        for (j = FLIP (head) ; j != EMPTY && Next[j] > EMPTY ; j = Next [j])
        {
        /* compare j with all nodes i following it in hash bucket */
        ASSERT (j >= 0 && j < n && Hash [j] == hash) ;
        p = Cp [j] ;
        pend = Cp [j+1] ;
        jlen = pend - p ;
        jscattered = FALSE ;
        DEBUG (for (i = 0 ; i < n ; i++) ASSERT (!Cmap_MARKED (i))) ;
        DEBUG (pruned = FALSE) ;
        ilast = j ;
        for (i = Next [j] ; i != EMPTY ; i = Next [i])
        {
            ASSERT (i >= 0 && i < n && Hash [i] == hash && i != j) ;
            pi = Cp [i] ;
            piend = Cp [i+1] ;
            ilen = piend - pi ;
            DEBUG (work++) ;
            if (ilen != jlen)
            {
            /* i and j have different degrees */
            ilast = i ;
            continue ;
            }
            /* scatter the pattern of node j, if not already */
            if (!jscattered)
            {
            Cmap_MARK (j) ;
            for ( ; p < pend ; p++)
            {
                Cmap_MARK (Ci [p]) ;
            }
            jscattered = TRUE ;
            DEBUG (work += jlen) ;
            }
            for (ok = Cmap_MARKED (i) ; ok && pi < piend ; pi++)
            {
            ok = Cmap_MARKED (Ci [pi]) ;
            DEBUG (work++) ;
            }
            if (ok)
            {
            /* found it.  kill node i and merge it into j */
            PRINT2 (("found "ID" absorbed into "ID"\n", i, j)) ;
            Hash [i] = FLIP (j) ;
            Cnw [j] += Cnw [i] ;
            Cnw [i] = 0 ;
            ASSERT (ilast != i && ilast >= 0 && ilast < n) ;
            Next [ilast] = Next [i] ; /* delete i from bucket */
            nodes_pruned++ ;
            DEBUG (goodwork += (ilen+1)) ;
            DEBUG (pruned = TRUE) ;
            }
            else
            {
            /* i and j are different */
            ilast = i ;
            }
        }
        DEBUG (if (pruned) goodwork += jlen) ;
        }
        /* empty the hash bucket, restoring Cew */
        Hhead [hash] = 1 ;
    }

    DEBUG (if (((work - goodwork) / (double) nz) > 0.20) PRINT0 ((
        "work %12g good %12g nz %12g (wasted work/nz: %6.2f )\n",
        work, goodwork, (double) nz, (work - goodwork) / ((double) nz)))) ;

    /* All hash buckets now empty.  Cmap no longer needed as workspace. ]
     * Cew no longer needed as Hhead; Cew is now restored to all ones. ]
     * Part no longer needed as workspace for Next. ] */
    }

    /* Edge weights are all one, node weights reflect node absorption */
    DEBUG (for (p = 0 ; p < csize ; p++) ASSERT (Cew [p] == 1)) ;
    DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) cnt += Cnw [j]) ;
    ASSERT (cnt == total_weight) ;

    /* ---------------------------------------------------------------------- */
    /* compress and partition the graph */
    /* ---------------------------------------------------------------------- */

    if (nodes_pruned == 0)
    {

    /* ------------------------------------------------------------------ */
    /* no pruning done at all.  Do not create the compressed graph */
    /* ------------------------------------------------------------------ */

    /* FUTURE WORK: could call CHACO, SCOTCH, ... here too */
    csep = CHOLMOD(metis_bisector) (C, Cnw, Cew, Part, Common) ;

    }
    else if (nodes_pruned == n-1)
    {

    /* ------------------------------------------------------------------ */
    /* only one node left.  This is a dense graph */
    /* ------------------------------------------------------------------ */

    PRINT2 (("completely dense graph\n")) ;
    csep = total_weight ;
    for (j = 0 ; j < n ; j++)
    {
        Part [j] = 2 ;
    }

    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* compress the graph and partition the compressed graph */
    /* ------------------------------------------------------------------ */

    /* ------------------------------------------------------------------ */
    /* create the map from the uncompressed graph to the compressed graph */
    /* ------------------------------------------------------------------ */

    /* Cmap [j] = k if node j is alive and the kth node of compressed graph.
     * The mapping is done monotonically (that is, k <= j) to simplify the
     * uncompression later on.  Cmap [j] = EMPTY if node j is dead. */

    for (j = 0 ; j < n ; j++)
    {
        Cmap [j] = EMPTY ;
    }
    k = 0 ;
    for (j = 0 ; j < n ; j++)
    {
        if (Cnw [j] > 0)
        {
        ASSERT (k <= j) ;
        Cmap [j] = k++ ;
        }
    }
    cn = k ;        /* # of nodes in compressed graph */
    PRINT2 (("compressed graph from "ID" to "ID" nodes\n", n, cn)) ;
    ASSERT (cn > 1 && cn == n - nodes_pruned) ;

    /* ------------------------------------------------------------------ */
    /* create the compressed graph */
    /* ------------------------------------------------------------------ */

    k = 0 ;
    pdest = 0 ;
    for (j = 0 ; j < n ; j++)
    {
        if (Cnw [j] > 0)
        {
        /* node j in the full graph is node k in the compressed graph */
        ASSERT (k <= j && Cmap [j] == k) ;
        p = Cp [j] ;
        pend = Cp [j+1] ;
        Cp [k] = pdest ;
        Cnw [k] = Cnw [j] ;
        for ( ; p < pend ; p++)
        {
            /* prune dead nodes, and remap to new node numbering */
            i = Ci [p] ;
            ASSERT (i >= 0 && i < n && i != j) ;
            i = Cmap [i] ;
            ASSERT (i >= EMPTY && i < cn && i != k) ;
            if (i > EMPTY)
            {
            ASSERT (pdest <= p) ;
            Ci [pdest++] = i ;
            }
        }
        k++ ;
        }
    }
    Cp [cn] = pdest ;
    C->nrow = cn ;
    C->ncol = cn ;  /* affects mem stats unless restored when C free'd */

#ifndef NDEBUG
    PRINT2 (("pruned graph ("ID"/"ID") nodes, ("ID"/"ID") edges\n",
            cn, n, pdest, nz)) ;
    PRINT2 (("compressed graph:\n")) ;
    for (cnt = 0, j = 0 ; j < cn ; j++)
    {
        PRINT2 ((""ID": ", j)) ;
        for (p = Cp [j] ; p < Cp [j+1] ; p++)
        {
        i = Ci [p] ;
        PRINT3 ((""ID" ", i)) ;
        ASSERT (i >= 0 && i < cn && i != j) ;
        }
        PRINT2 (("weight: "ID"\n", Cnw [j])) ;
        ASSERT (Cnw [j] > 0) ;
        cnt += Cnw [j] ;
    }
    ASSERT (cnt == total_weight) ;
    for (j = 0 ; j < n ; j++) PRINT2 (("Cmap ["ID"] = "ID"\n", j, Cmap[j]));
    ASSERT (k == cn) ;
#endif

    /* ------------------------------------------------------------------ */
    /* find the separator of the compressed graph */
    /* ------------------------------------------------------------------ */

    /* FUTURE WORK: could call CHACO, SCOTCH, ... here too */
    csep = CHOLMOD(metis_bisector) (C, Cnw, Cew, Part, Common) ;

    if (csep < 0)
    {
        /* failed */
        return (-1) ;
    }

    PRINT2 (("Part: ")) ;
    DEBUG (for (j = 0 ; j < cn ; j++) PRINT2 ((""ID" ", Part [j]))) ;
    PRINT2 (("\n")) ;

    /* Cp and Ci no longer needed */

    /* ------------------------------------------------------------------ */
    /* find the separator of the uncompressed graph */
    /* ------------------------------------------------------------------ */

    /* expand the separator to live nodes in the uncompressed graph */
    for (j = n-1 ; j >= 0 ; j--)
    {
        /* do this in reverse order so that Cnw can be expanded in place */
        k = Cmap [j] ;
        ASSERT (k >= EMPTY && k < n) ;
        if (k > EMPTY)
        {
        /* node k in compressed graph and is node j in full graph */
        ASSERT (k <= j) ;
        ASSERT (Hash [j] >= EMPTY) ;
        Part [j] = Part [k] ;
        Cnw [j] = Cnw [k] ;
        }
        else
        {
        /* node j is a dead node */
        Cnw [j] = 0 ;
        DEBUG (Part [j] = EMPTY) ;
        ASSERT (Hash [j] < EMPTY) ;
        }
    }

    /* find the components for the dead nodes */
    for (i = 0 ; i < n ; i++)
    {
        if (Hash [i] < EMPTY)
        {
        /* node i has been absorbed into node j */
        j = FLIP (Hash [i]) ;
        ASSERT (Part [i] == EMPTY && j >= 0 && j < n && Cnw [i] == 0) ;
        Part [i] = Part [j] ;
        }
        ASSERT (Part [i] >= 0 && Part [i] <= 2) ;
    }

#ifndef NDEBUG
    PRINT2 (("Part: ")) ;
    for (cnt = 0, j = 0 ; j < n ; j++)
    {
        ASSERT (Part [j] != EMPTY) ;
        PRINT2 ((""ID" ", Part [j])) ;
        if (Part [j] == 2) cnt += Cnw [j] ;
    }
    PRINT2 (("\n")) ;
    PRINT2 (("csep "ID" "ID"\n", cnt, csep)) ;
    ASSERT (cnt == csep) ;
    for (cnt = 0, j = 0 ; j < n ; j++) cnt += Cnw [j] ;
    ASSERT (cnt == total_weight) ;
#endif

    }

    /* ---------------------------------------------------------------------- */
    /* return the separator (or -1 if error) */
    /* ---------------------------------------------------------------------- */

    PRINT2 (("Partition done, n "ID" csep "ID"\n", n, csep)) ;
    return (csep) ;
}


/* ========================================================================== */
/* === clear_flag =========================================================== */
/* ========================================================================== */

/* A node j has been removed from the graph if Flag [j] < EMPTY.
 * If Flag [j] >= EMPTY && Flag [j] < mark, then node j is alive but unmarked.
 * Flag [j] == mark means that node j is alive and marked.  Incrementing mark
 * means that all nodes are either (still) dead, or live but unmarked.
 *
 * If Map is NULL, then on output, Common->mark < Common->Flag [i] for all i
 * from 0 to Common->nrow.  This is the same output condition as
 * cholmod_clear_flag, except that this routine maintains the Flag [i] < EMPTY
 * condition as well, if that condition was true on input.
 *
 * If Map is non-NULL, then on output, Common->mark < Common->Flag [i] for all
 * i in the set Map [0..cn-1].
 *
 * workspace: Flag (nrow)
 */

static UF_long clear_flag (Int *Map, Int cn, cholmod_common *Common)
{
    Int nrow, i ;
    Int *Flag ;
    PRINT2 (("old mark %ld\n", Common->mark)) ;
    Common->mark++ ;
    PRINT2 (("new mark %ld\n", Common->mark)) ;
    if (Common->mark <= 0)
    {
    nrow = Common->nrow ;
    Flag = Common->Flag ;
        if (Map != NULL)
        {
            for (i = 0 ; i < cn ; i++)
            {
                /* if Flag [Map [i]] < EMPTY, leave it alone */
                if (Flag [Map [i]] >= EMPTY)
                {
                    Flag [Map [i]] = EMPTY ;
                }
            }
            /* now Flag [Map [i]] <= EMPTY for all i */
        }
        else
        {
            for (i = 0 ; i < nrow ; i++)
            {
                /* if Flag [i] < EMPTY, leave it alone */
                if (Flag [i] >= EMPTY)
                {
                    Flag [i] = EMPTY ;
                }
            }
            /* now Flag [i] <= EMPTY for all i */
        }
    Common->mark = 0 ;
    }
    return (Common->mark) ;
}


/* ========================================================================== */
/* === find_components ====================================================== */
/* ========================================================================== */

/* Find all connected components of the current subgraph C.  The subgraph C
 * consists of the nodes of B that appear in the set Map [0..cn-1].  If Map
 * is NULL, then it is assumed to be the identity mapping
 * (Map [0..cn-1] = 0..cn-1).
 *
 * A node j does not appear in B if it has been ordered (Flag [j] < EMPTY,
 * which means that j has been ordered and is "deleted" from B).
 *
 * If the size of a component is large, it is placed on the component stack,
 * Cstack.  Otherwise, its nodes are ordered and it is not placed on the Cstack.
 *
 * A component S is defined by a "representative node" (repnode for short)
 * called the snode, which is one of the nodes in the subgraph.  Likewise, the
 * subgraph C is defined by its repnode, called cnode.
 * 
 * If Part is not NULL on input, then Part [i] determines how the components
 * are placed on the stack.  Components containing nodes i with Part [i] == 0
 * are placed first, followed by components with Part [i] == 1. 
 *
 * The first node placed in each of the two parts is flipped when placed in
 * the Cstack.  This allows the components of the two parts to be found simply
 * by traversing the Cstack.
 *
 * workspace: Flag (nrow)
 */

static void find_components
(
    /* inputs, not modified on output */
    cholmod_sparse *B,
    Int Map [ ],        /* size n, only Map [0..cn-1] used */
    Int cn,         /* # of nodes in C */
    Int cnode,          /* root node of component C, or EMPTY if C is the
                 * entire graph B */

    Int Part [ ],       /* size cn, optional */

    /* input/output */
    Int Bnz [ ],        /* size n.  Bnz [j] = # nonzeros in column j of B.
                 * Reduce since B is pruned of dead nodes. */

    Int CParent [ ],        /* CParent [i] = j if component with repnode j is
                 * the parent of the component with repnode i.
                 * CParent [i] = EMPTY if the component with
                 * repnode i is a root of the separator tree.
                 * CParent [i] is -2 if i is not a repnode. */
    Int Cstack [ ],     /* component stack for nested dissection */
    Int *top,           /* Cstack [0..top] contains root nodes of the
                 * the components currently in the stack */

    /* workspace, undefined on input and output: */
    Int Queue [ ],      /* size n, for breadth-first search */

    cholmod_common *Common
)
{
    Int n, mark, cj, j, sj, sn, p, i, snode, pstart, pdest, pend, nd_components,
    part, first, save_mark ;
    Int *Bp, *Bi, *Flag ;

    /* ---------------------------------------------------------------------- */
    /* get workspace */
    /* ---------------------------------------------------------------------- */

    PRINT2 (("find components: cn %d\n", cn)) ;
    Flag = Common->Flag ;       /* size n */

    /* force initialization of Flag [Map [0..cn-1]] */
    save_mark = Common->mark ;      /* save the current mark */
    Common->mark = EMPTY ;

    /* clear Flag; preserve Flag [Map [i]] if Flag [Map [i]] already < EMPTY */
    /* this takes O(cn) time */
    mark = clear_flag (Map, cn, Common) ;

    Bp = B->p ;
    Bi = B->i ;
    n = B->nrow ;
    ASSERT (cnode >= EMPTY && cnode < n) ;
    ASSERT (IMPLIES (cnode >= 0, Flag [cnode] < EMPTY)) ;

    /* get ordering parameters */
    nd_components = Common->method [Common->current].nd_components ;

    /* ---------------------------------------------------------------------- */
    /* find the connected components of C via a breadth-first search */
    /* ---------------------------------------------------------------------- */

    part = (Part == NULL) ? 0 : 1 ;

    /* examine each part (part 1 and then part 0) */
    for (part = (Part == NULL) ? 0 : 1 ; part >= 0 ; part--)
    {

    /* first is TRUE for the first connected component in each part */
    first = TRUE ;

    /* find all connected components in the current part */
    for (cj = 0 ; cj < cn ; cj++)
    {
        /* get node snode, which is node cj of C.  It might already be in
         * the separator of C (and thus ordered, with Flag [snode] < EMPTY)
         */
        snode = (Map == NULL) ? (cj) : (Map [cj]) ;
        ASSERT (snode >= 0 && snode < n) ;

        if (Flag [snode] >= EMPTY && Flag [snode] < mark
            && ((Part == NULL) || Part [cj] == part))
        {

        /* ---------------------------------------------------------- */
        /* find new connected component S */
        /* ---------------------------------------------------------- */

        /* node snode is the repnode of a connected component S, the
         * parent of which is cnode, the repnode of C.  If cnode is
         * EMPTY then C is the original graph B. */
        PRINT2 (("----------:::snode "ID" cnode "ID"\n", snode, cnode));

        ASSERT (CParent [snode] == -2) ;
        if (first || nd_components)
        {
            /* If this is the first node in this part, then it becomes
             * the repnode of all components in this part, and all
             * components in this part form a single node in the
             * separator tree.  If nd_components is TRUE, then all
             * connected components form their own node in the
             * separator tree.
             */
            CParent [snode] = cnode ;
        }

        /* place j in the queue and mark it */
        Queue [0] = snode ;
        Flag [snode] = mark ;
        sn = 1 ;

        /* breadth-first traversal, starting at node j */
        for (sj = 0 ; sj < sn ; sj++)
        {
            /* get node j from head of Queue and traverse its edges */
            j = Queue [sj] ;
            PRINT2 (("    j: "ID"\n", j)) ;
            ASSERT (j >= 0 && j < n) ;
            ASSERT (Flag [j] == mark) ;
            pstart = Bp [j] ;
            pdest = pstart ;
            pend = pstart + Bnz [j] ;
            for (p = pstart ; p < pend ; p++)
            {
            i = Bi [p] ;
            if (i != j && Flag [i] >= EMPTY)
            {
                /* node is still in the graph */
                Bi [pdest++] = i ;
                if (Flag [i] < mark)
                {
                /* node i is in this component S, and unflagged
                 * (first time node i has been seen in this BFS)
                 * place node i in the queue and mark it */
                Queue [sn++] = i ;
                Flag [i] = mark ;
                }
            }
            }
            /* edges to dead nodes have been removed */
            Bnz [j] = pdest - pstart ;
        }

        /* ---------------------------------------------------------- */
        /* order S if it is small; place it on Cstack otherwise */
        /* ---------------------------------------------------------- */

        PRINT2 (("sn "ID"\n", sn)) ;

        /* place the new component on the Cstack.  Flip the node if
         * is the first connected component of the current part,
         * or if all components are treated as their own node in
         * the separator tree. */
        Cstack [++(*top)] =
            (first || nd_components) ? FLIP (snode) : snode ;
        first = FALSE ;
        }
    }
    }

    /* restore the flag (normally taking O(1) time except for Int overflow) */
    Common->mark = save_mark++ ;
    clear_flag (NULL, 0, Common) ;
    DEBUG (for (i = 0 ; i < n ; i++) ASSERT (Flag [i] < Common->mark)) ;
}


/* ========================================================================== */
/* === cholmod_bisect ======================================================= */
/* ========================================================================== */

/* Finds a node bisector of A, A*A', A(:,f)*A(:,f)'.
 *
 * workspace: Flag (nrow),
 *  Iwork (nrow if symmetric, max (nrow,ncol) if unsymmetric).
 *  Allocates a temporary matrix B=A*A' or B=A,
 *  and O(nnz(A)) temporary memory space.
 */

UF_long CHOLMOD(bisect) /* returns # of nodes in separator */
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to bisect */
    Int *fset,      /* subset of 0:(A->ncol)-1 */
    size_t fsize,   /* size of fset */
    int compress,   /* if TRUE, compress the graph first */
    /* ---- output --- */
    Int *Partition, /* size A->nrow.  Node i is in the left graph if
             * Partition [i] = 0, the right graph if 1, and in the
             * separator if 2. */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Bp, *Bi, *Hash, *Cmap, *Bnw, *Bew, *Iwork ;
    cholmod_sparse *B ;
    unsigned Int hash ;
    Int j, n, bnz, sepsize, p, pend ;
    size_t csize, s ;
    int ok = TRUE ;

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (EMPTY) ;
    RETURN_IF_NULL (A, EMPTY) ;
    RETURN_IF_NULL (Partition, EMPTY) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
    Common->status = CHOLMOD_OK ;

    /* ---------------------------------------------------------------------- */
    /* quick return */
    /* ---------------------------------------------------------------------- */

    n = A->nrow ;
    if (n == 0)
    {
    return (0) ;
    }

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    /* s = n + MAX (n, A->ncol) */
    s = CHOLMOD(add_size_t) (A->nrow, MAX (A->nrow, A->ncol), &ok) ;
    if (!ok)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    return (EMPTY) ;
    }

    CHOLMOD(allocate_work) (n, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (EMPTY) ;
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

    Iwork = Common->Iwork ;
    Hash = Iwork ;      /* size n, (i/l/l) */
    Cmap = Iwork + n ;      /* size n, (i/i/l) */

    /* ---------------------------------------------------------------------- */
    /* convert the matrix to adjacency list form */
    /* ---------------------------------------------------------------------- */

    /* The input graph to must be symmetric, with no diagonal entries
     * present.  The columns need not be sorted. */

    /* B = A, A*A', or A(:,f)*A(:,f)', upper and lower parts present */

    if (A->stype)
    {
    /* Add the upper/lower part to a symmetric lower/upper matrix by
     * converting to unsymmetric mode */
    /* workspace: Iwork (nrow) */
    B = CHOLMOD(copy) (A, 0, -1, Common) ;
    }
    else
    {
    /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
    /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
    B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
    }

    if (Common->status < CHOLMOD_OK)
    {
    return (EMPTY) ;
    }
    Bp = B->p ;
    Bi = B->i ;
    bnz = Bp [n] ;
    ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;

    /* B does not include the diagonal, and both upper and lower parts.
     * Common->anz includes the diagonal, and just the lower part of B */
    Common->anz = bnz / 2 + ((double) n) ;

    /* Bew should be at least size n for the hash function to work well */
    /* this cannot cause overflow, because the matrix is already created */
    csize = MAX (((size_t) n) + 1, (size_t) bnz) ;

    /* create the graph using Flag as workspace for node weights [ */
    Bnw = Common->Flag ;    /* size n workspace */

    /* compute hash for each node if compression requested */
    if (compress)
    {
    for (j = 0 ; j < n ; j++)
    {
        hash = j ;
        pend = Bp [j+1] ;
        for (p = Bp [j] ; p < pend ; p++)
        {
        hash += Bi [p] ;
        ASSERT (Bi [p] != j) ;
        }
        /* finalize the hash key for node j */
        hash %= csize ;
        Hash [j] = (Int) hash ;
        ASSERT (Hash [j] >= 0 && Hash [j] < csize) ;
    }
    }

    /* allocate edge weights */
    Bew = CHOLMOD(malloc) (csize, sizeof (Int), Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    /* out of memory */
    CHOLMOD(free_sparse) (&B, Common) ;
    CHOLMOD(free) (csize, sizeof (Int), Bew, Common) ;
    return (EMPTY) ;
    }

    /* graph has unit node and edge weights */
    for (j = 0 ; j < n ; j++)
    {
    Bnw [j] = 1 ;
    }
    for (s = 0 ; s < csize ; s++)
    {
    Bew [s] = 1 ;
    }

    /* ---------------------------------------------------------------------- */
    /* compress and partition the graph */
    /* ---------------------------------------------------------------------- */

    sepsize = partition (
#ifndef NDEBUG
        csize,
#endif
        compress, Hash, B, Bnw, Bew, Cmap, Partition, Common) ;

    /* contents of Bp, Bi, Bnw, and Bew no longer needed ] */

    /* If partition fails, free the workspace below and return sepsize < 0 */

    /* ---------------------------------------------------------------------- */
    /* free workspace */
    /* ---------------------------------------------------------------------- */

    B->ncol = n ;   /* restore size for memory usage statistics */
    CHOLMOD(free_sparse) (&B, Common) ;
    Common->mark = EMPTY ;
    CHOLMOD_CLEAR_FLAG (Common) ;
    CHOLMOD(free) (csize, sizeof (Int), Bew, Common) ;
    return (sepsize) ;
}


/* ========================================================================== */
/* === cholmod_nested_dissection ============================================ */
/* ========================================================================== */

/* This method uses a node bisector, applied recursively (but using a
 * non-recursive algorithm).  Once the graph is partitioned, it calls a
 * constrained min degree code (CAMD or CSYMAMD for A+A', and CCOLAMD for A*A')
 * to order all the nodes in the graph - but obeying the constraints determined
 * by the separators.  This routine is similar to METIS_NodeND, except for how
 * it treats the leaf nodes.  METIS_NodeND orders the leaves of the separator
 * tree with MMD, ignoring the rest of the matrix when ordering a single leaf.
 * This routine orders the whole matrix with CSYMAMD or CCOLAMD, all at once,
 * when the graph partitioning is done.
 *
 * This function also returns a postorderd separator tree (CParent), and a
 * mapping of nodes in the graph to nodes in the separator tree (Cmember).
 *
 * workspace: Flag (nrow), Head (nrow+1), Iwork (4*nrow + (ncol if unsymmetric))
 *  Allocates a temporary matrix B=A*A' or B=A,
 *  and O(nnz(A)) temporary memory space.
 *  Allocates an additional 3*n*sizeof(Int) temporary workspace
 */

UF_long CHOLMOD(nested_dissection) /* returns # of components, or -1 if error */
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to order */
    Int *fset,      /* subset of 0:(A->ncol)-1 */
    size_t fsize,   /* size of fset */
    /* ---- output --- */
    Int *Perm,      /* size A->nrow, output permutation */
    Int *CParent,   /* size A->nrow.  On output, CParent [c] is the parent
             * of component c, or EMPTY if c is a root, and where
             * c is in the range 0 to # of components minus 1 */
    Int *Cmember,   /* size A->nrow.  Cmember [j] = c if node j of A is
             * in component c */
    /* --------------- */
    cholmod_common *Common
)
{
    double prune_dense, nd_oksep ;
    Int *Bp, *Bi, *Bnz, *Cstack, *Imap, *Map, *Flag, *Head, *Next, *Bnw, *Iwork,
    *Ipost, *NewParent, *Hash, *Cmap, *Cp, *Ci, *Cew, *Cnw, *Part, *Post,
    *Work3n ;
    unsigned Int hash ;
    Int n, bnz, top, i, j, k, cnode, cdense, p, cj, cn, ci, cnz, mark, c, uncol,
    sepsize, parent, ncomponents, threshold, ndense, pstart, pdest, pend,
    nd_compress, nd_camd, csize, jnext, nd_small, total_weight,
    nchild, child = EMPTY ;
    cholmod_sparse *B, *C ;
    size_t s ;
    int ok = TRUE ;
    DEBUG (Int cnt) ;

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (EMPTY) ;
    RETURN_IF_NULL (A, EMPTY) ;
    RETURN_IF_NULL (Perm, EMPTY) ;
    RETURN_IF_NULL (CParent, EMPTY) ;
    RETURN_IF_NULL (Cmember, EMPTY) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
    Common->status = CHOLMOD_OK ;

    /* ---------------------------------------------------------------------- */
    /* quick return */
    /* ---------------------------------------------------------------------- */

    n = A->nrow ;
    if (n == 0)
    {
    return (1) ;
    }

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    /* get ordering parameters */
    prune_dense = Common->method [Common->current].prune_dense ;
    nd_compress = Common->method [Common->current].nd_compress ;
    nd_oksep = Common->method [Common->current].nd_oksep ;
    nd_oksep = MAX (0, nd_oksep) ;
    nd_oksep = MIN (1, nd_oksep) ;
    nd_camd = Common->method [Common->current].nd_camd ;
    nd_small = Common->method [Common->current].nd_small ;
    nd_small = MAX (4, nd_small) ;

    PRINT0 (("nd_components %d nd_small %d nd_oksep %g\n", 
    Common->method [Common->current].nd_components,
    nd_small, nd_oksep)) ;

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    /* s = 4*n + uncol */
    uncol = (A->stype == 0) ? A->ncol : 0 ;
    s = CHOLMOD(mult_size_t) (n, 4, &ok) ;
    s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
    if (!ok)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    return (EMPTY) ;
    }

    CHOLMOD(allocate_work) (n, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (EMPTY) ;
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

    /* ---------------------------------------------------------------------- */
    /* get workspace */
    /* ---------------------------------------------------------------------- */

    Flag = Common->Flag ;   /* size n */
    Head = Common->Head ;   /* size n+1, all equal to -1 */

    Iwork = Common->Iwork ;
    Imap = Iwork ;      /* size n, same as Queue in find_components */
    Map  = Iwork + n ;      /* size n */
    Bnz  = Iwork + 2*((size_t) n) ; /* size n */
    Hash = Iwork + 3*((size_t) n) ; /* size n */

    Work3n = CHOLMOD(malloc) (n, 3*sizeof (Int), Common) ;
    Part = Work3n ;     /* size n */
    Bnw  = Part + n ;       /* size n */
    Cnw  = Bnw + n ;        /* size n */

    Cstack = Perm ;     /* size n, use Perm as workspace for Cstack [ */
    Cmap = Cmember ;        /* size n, use Cmember as workspace [ */

    if (Common->status < CHOLMOD_OK)
    {
    return (EMPTY) ;
    }

    /* ---------------------------------------------------------------------- */
    /* convert B to symmetric form with both upper/lower parts present */
    /* ---------------------------------------------------------------------- */

    /* B = A+A', A*A', or A(:,f)*A(:,f)', upper and lower parts present */

    if (A->stype)
    {
    /* Add the upper/lower part to a symmetric lower/upper matrix by
     * converting to unsymmetric mode */
    /* workspace: Iwork (nrow) */
    B = CHOLMOD(copy) (A, 0, -1, Common) ;
    }
    else
    {
    /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
    /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
    B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
    }

    if (Common->status < CHOLMOD_OK)
    {
    CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
    return (EMPTY) ;
    }
    Bp = B->p ;
    Bi = B->i ;
    bnz = CHOLMOD(nnz) (B, Common) ;
    ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;
    csize = MAX (n, bnz) ;
    ASSERT (CHOLMOD(dump_sparse) (B, "B for nd:", Common) >= 0) ;

    /* ---------------------------------------------------------------------- */
    /* initializations */
    /* ---------------------------------------------------------------------- */

    /* all nodes start out unmarked and unordered (Type 4, see below) */
    Common->mark = EMPTY ;
    CHOLMOD_CLEAR_FLAG (Common) ;
    ASSERT (Flag == Common->Flag) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

    for (j = 0 ; j < n ; j++)
    {
    CParent [j] = -2 ;
    }

    /* prune dense nodes from B */
    if (IS_NAN (prune_dense) || prune_dense < 0)
    {
    /* only remove completely dense nodes */
    threshold = n-2 ;
    }
    else
    {
    /* remove nodes with degree more than threshold */
    threshold = (Int) (MAX (16, prune_dense * sqrt ((double) (n)))) ;
    threshold = MIN (n, threshold) ;
    }
    ndense = 0 ;
    cnode = EMPTY ;
    cdense = EMPTY ;

    for (j = 0 ; j < n ; j++)
    {
    Bnz [j] = Bp [j+1] - Bp [j] ;
    if (Bnz [j] > threshold)
    {
        /* node j is dense, prune it from B */
        PRINT2 (("j is dense %d\n", j)) ;
        ndense++ ;
        if (cnode == EMPTY)
        {
        /* first dense node found becomes root of this component,
         * which contains all of the dense nodes found here */
        cdense = j ;
        cnode = j ;
        CParent [cnode] = EMPTY ;
        }
        Flag [j] = FLIP (cnode) ;
    }
    }
    B->packed = FALSE ;
    ASSERT (B->nz == NULL) ;

    if (ndense == n)
    {
    /* all nodes removed: Perm is identity, all nodes in component zero,
     * and the separator tree has just one node. */
    PRINT2 (("all nodes are dense\n")) ;
    for (k = 0 ; k < n ; k++)
    {
        Perm [k] = k ;
        Cmember [k] = 0 ;
    }
    CParent [0] = EMPTY ;
    CHOLMOD(free_sparse) (&B, Common) ;
    CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
    Common->mark = EMPTY ;
    CHOLMOD_CLEAR_FLAG (Common) ;
    return (1) ;
    }

    /* Cp and Ci are workspace to construct the subgraphs to partition */
    C = CHOLMOD(allocate_sparse) (n, n, csize, FALSE, TRUE, 0, CHOLMOD_PATTERN,
        Common) ;
    Cew  = CHOLMOD(malloc) (csize, sizeof (Int), Common) ;

    if (Common->status < CHOLMOD_OK)
    {
    /* out of memory */
    CHOLMOD(free_sparse) (&C, Common) ;
    CHOLMOD(free_sparse) (&B, Common) ;
    CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
    CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
    Common->mark = EMPTY ;
    CHOLMOD_CLEAR_FLAG (Common) ;
    PRINT2 (("out of memory for C, etc\n")) ;
    return (EMPTY) ;
    }

    Cp = C->p ;
    Ci = C->i ;

    /* create initial unit node and edge weights */
    for (j = 0 ; j < n ; j++)
    {
    Bnw [j] = 1 ;
    }
    for (p = 0 ; p < csize ; p++)
    {
    Cew [p] = 1 ;
    }

    /* push the initial connnected components of B onto the Cstack */
    top = EMPTY ;   /* Cstack is empty */
    /* workspace: Flag (nrow), Iwork (nrow); use Imap as workspace for Queue [*/
    find_components (B, NULL, n, cnode, NULL,
        Bnz, CParent, Cstack, &top, Imap, Common) ;
    /* done using Imap as workspace for Queue ] */

    /* Nodes can now be of Type 0, 1, 2, or 4 (see definition below) */

    /* ---------------------------------------------------------------------- */
    /* while Cstack is not empty, do: */
    /* ---------------------------------------------------------------------- */

    while (top >= 0)
    {

    /* clear the Flag array, but do not modify negative entries in Flag  */
    mark = clear_flag (NULL, 0, Common) ;

    DEBUG (for (i = 0 ; i < n ; i++) Imap [i] = EMPTY) ;

    /* ------------------------------------------------------------------ */
    /* get node(s) from the top of the Cstack */
    /* ------------------------------------------------------------------ */

    /* i is the repnode of its (unordered) connected component.  Get
     * all repnodes for all connected components of a single part.  If
     * each connected component is to be ordered separately (nd_components
     * is TRUE), then this while loop iterates just once. */

    cnode = EMPTY ;
    cn = 0 ;
    while (cnode == EMPTY)
    {
        i = Cstack [top--] ;

        if (i < 0)
        {
        /* this is the last node in this component */
        i = FLIP (i) ;
        cnode = i ;
        }

        ASSERT (i >= 0 && i < n && Flag [i] >= EMPTY) ;

        /* place i in the queue and mark it */
        Map [cn] = i ;
        Flag [i] = mark ;
        Imap [i] = cn ;
        cn++ ;
    }

    ASSERT (cnode != EMPTY) ;

    /* During ordering, there are five kinds of nodes in the graph of B,
     * based on Flag [j] and CParent [j] for nodes j = 0 to n-1:
     *
     * Type 0: If cnode is a repnode of an unordered component, then
     * CParent [cnode] is in the range EMPTY to n-1 and
     * Flag [cnode] >= EMPTY.  This is a "live" node.
     *
     * Type 1: If cnode is a repnode of an ordered separator component,
     * then Flag [cnode] < EMPTY and FLAG [cnode] = FLIP (cnode).
     * CParent [cnode] is in the range EMPTY to n-1.  cnode is a root of
     * the separator tree if CParent [cnode] == EMPTY.  This node is dead.
     *
     * Type 2: If node j isn't a repnode, has not been absorbed via
     * graph compression into another node, but is in an ordered separator
     * component, then cnode = FLIP (Flag [j]) gives the repnode of the
     * component that contains j and CParent [j]  is -2.  This node is dead.
     * Note that Flag [j] < EMPTY.
     *
     * Type 3: If node i has been absorbed via graph compression into some
     * other node j = FLIP (Flag [i]) where j is not a repnode.
     * CParent [j] is -2.  Node i may or may not be in an ordered
     * component.  This node is dead.  Note that Flag [j] < EMPTY.
     *
     * Type 4: If node j is "live" (not in an ordered component, and not
     * absorbed into any other node), then Flag [j] >= EMPTY.
     *
     * Only "live" nodes (of type 0 or 4) are placed in a subgraph to be
     * partitioned.  Node j is alive if Flag [j] >= EMPTY, and dead if
     * Flag [j] < EMPTY.
     */

    /* ------------------------------------------------------------------ */
    /* create the subgraph for this connected component C */
    /* ------------------------------------------------------------------ */

    /* Do a breadth-first search of the graph starting at cnode.
     * use Map [0..cn-1] for nodes in the component C [
     * use Cnw and Cew for node and edge weights of the resulting subgraph [
     * use Cp and Ci for the resulting subgraph [
     * use Imap [i] for all nodes i in B that are in the component C [
     */

    cnz = 0 ;
    total_weight = 0 ;
    for (cj = 0 ; cj < cn ; cj++)
    {
        /* get node j from the head of the queue; it is node cj of C */
        j = Map [cj] ;
        ASSERT (Flag [j] == mark) ;
        Cp [cj] = cnz ;
        Cnw [cj] = Bnw [j] ;
        ASSERT (Cnw [cj] >= 0) ;
        total_weight += Cnw [cj] ;
        pstart = Bp [j] ;
        pdest = pstart ;
        pend = pstart + Bnz [j] ;
        hash = cj ;
        for (p = pstart ; p < pend ; p++)
        {
        i = Bi [p] ;
        /* prune diagonal entries and dead edges from B */
        if (i != j && Flag [i] >= EMPTY)
        {
            /* live node i is in the current component */
            Bi [pdest++] = i ;
            if (Flag [i] != mark)
            {
            /* First time node i has been seen, it is a new node
             * of C.  place node i in the queue and mark it */
            Map [cn] = i ;
            Flag [i] = mark ;
            Imap [i] = cn ;
            cn++ ;
            }
            /* place the edge (cj,ci) in the adjacency list of cj */
            ci = Imap [i] ;
            ASSERT (ci >= 0 && ci < cn && ci != cj && cnz < csize) ;
            Ci [cnz++] = ci ;
            hash += ci ;
        }
        }
        /* edges to dead nodes have been removed */
        Bnz [j] = pdest - pstart ;
        /* finalize the hash key for column j */
        hash %= csize ;
        Hash [cj] = (Int) hash ;
        ASSERT (Hash [cj] >= 0 && Hash [cj] < csize) ;
    }
    Cp [cn] = cnz ;
    C->nrow = cn ;
    C->ncol = cn ;  /* affects mem stats unless restored when C free'd */

    /* contents of Imap no longer needed ] */

#ifndef NDEBUG
    for (cj = 0 ; cj < cn ; cj++)
    {
        j = Map [cj] ;
        PRINT2 (("----------------------------C column cj: "ID" j: "ID"\n",
        cj, j)) ;
        ASSERT (j >= 0 && j < n) ;
        ASSERT (Flag [j] >= EMPTY) ;
        for (p = Cp [cj] ; p < Cp [cj+1] ; p++)
        {
        ci = Ci [p] ;
        i = Map [ci] ;
        PRINT3 (("ci: "ID" i: "ID"\n", ci, i)) ;
        ASSERT (ci != cj && ci >= 0 && ci < cn) ;
        ASSERT (i != j && i >= 0 && i < n) ;
        ASSERT (Flag [i] >= EMPTY) ;
        }
    }
#endif

    PRINT0 (("consider cn %d nd_small %d ", cn, nd_small)) ;
    if (cn < nd_small)  /* could be 'total_weight < nd_small' instead */
    {
        /* place all nodes in the separator */
        PRINT0 ((" too small\n")) ;
        sepsize = total_weight ;
    }
    else
    {

        /* Cp and Ci now contain the component, with cn nodes and cnz
         * nonzeros.  The mapping of a node cj into node j the main graph
         * B is given by Map [cj] = j */
        PRINT0 ((" cut\n")) ;

        /* -------------------------------------------------------------- */
        /* compress and partition the graph C */
        /* -------------------------------------------------------------- */

        /* The edge weights Cew [0..csize-1] are all 1's on input to and
         * output from the partition routine. */

        sepsize = partition (
#ifndef NDEBUG
            csize,
#endif
            nd_compress, Hash, C, Cnw, Cew,
            Cmap, Part, Common) ;

        /* contents of Cp and Ci no longer needed ] */

        if (sepsize < 0)
        {
        /* failed */
        C->ncol = n ;   /* restore size for memory usage statistics */
        CHOLMOD(free_sparse) (&C, Common) ;
        CHOLMOD(free_sparse) (&B, Common) ;
        CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
        CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
        Common->mark = EMPTY ;
        CHOLMOD_CLEAR_FLAG (Common) ;
        return (EMPTY) ;
        }

        /* -------------------------------------------------------------- */
        /* compress B based on how C was compressed */
        /* -------------------------------------------------------------- */

        for (ci = 0 ; ci < cn ; ci++)
        {
        if (Hash [ci] < EMPTY)
        {
            /* ci is dead in C, having been absorbed into cj */
            cj = FLIP (Hash [ci]) ;
            PRINT2 (("In C, "ID" absorbed into "ID" (wgt now "ID")\n",
                ci, cj, Cnw [cj])) ;
            /* i is dead in B, having been absorbed into j */
            i = Map [ci] ;
            j = Map [cj] ;
            PRINT2 (("In B, "ID" (wgt "ID") => "ID" (wgt "ID")\n",
                i, Bnw [i], j, Bnw [j], Cnw [cj])) ;
            /* more than one node may be absorbed into j.  This is
             * accounted for in Cnw [cj].  Assign it here rather
             * than += Bnw [i] */
            Bnw [i] = 0 ;
            Bnw [j] = Cnw [cj] ;
            Flag [i] = FLIP (j) ;
        }
        }

        DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) cnt += Bnw [j]) ;
        ASSERT (cnt == n) ;
    }

    /* contents of Cnw [0..cn-1] no longer needed ] */

    /* ------------------------------------------------------------------ */
    /* order the separator, and stack the components when C is split */
    /* ------------------------------------------------------------------ */

    /* one more component has been found: either the separator of C,
     * or all of C */

    ASSERT (sepsize >= 0 && sepsize <= total_weight) ;

    PRINT0 (("sepsize %d tot %d : %8.4f ", sepsize, total_weight,
        ((double) sepsize) / ((double) total_weight))) ;

    if (sepsize == total_weight || sepsize == 0 ||
        sepsize > nd_oksep * total_weight)
    {
        /* Order the nodes in the component.  The separator is too large,
         * or empty.  Note that the partition routine cannot return a
         * sepsize of zero, but it can return a separator consisting of the
         * whole graph.  The "sepsize == 0" test is kept, above, in case the
         * partition routine changes.  In either case, this component
         * remains unsplit, and becomes a leaf of the separator tree. */
        PRINT2 (("cnode %d sepsize zero or all of graph: "ID"\n",
        cnode, sepsize)) ;
        for (cj = 0 ; cj < cn ; cj++)
        {
        j = Map [cj] ;
        Flag [j] = FLIP (cnode) ;
        PRINT2 (("      node cj: "ID" j: "ID" ordered\n", cj, j)) ;
        }
        ASSERT (Flag [cnode] == FLIP (cnode)) ;
        ASSERT (cnode != EMPTY && Flag [cnode] < EMPTY) ;
        PRINT0 (("discarded\n")) ;

    }
    else
    {

        /* Order the nodes in the separator of C and find a new repnode
         * cnode that is in the separator of C.  This requires the separator
         * to be non-empty. */
        PRINT0 (("sepsize not tiny: "ID"\n", sepsize)) ;
        parent = CParent [cnode] ;
        ASSERT (parent >= EMPTY && parent < n) ;
        CParent [cnode] = -2 ;
        cnode = EMPTY ;
        for (cj = 0 ; cj < cn ; cj++)
        {
        j = Map [cj] ;
        if (Part [cj] == 2)
        {
            /* All nodes in the separator become part of a component
             * whose repnode is cnode */
            PRINT2 (("node cj: "ID" j: "ID" ordered\n", cj, j)) ;
            if (cnode == EMPTY)
            {
            PRINT2(("------------new cnode: cj "ID" j "ID"\n",
                    cj, j)) ;
            cnode = j ;
            }
            Flag [j] = FLIP (cnode) ;
        }
        else
        {
            PRINT2 (("      node cj: "ID" j: "ID" not ordered\n",
                cj, j)) ;
        }
        }
        ASSERT (cnode != EMPTY && Flag [cnode] < EMPTY) ;
        ASSERT (CParent [cnode] == -2) ;
        CParent [cnode] = parent ;

        /* find the connected components when C is split, and push
         * them on the Cstack.  Use Imap as workspace for Queue. [ */
        /* workspace: Flag (nrow) */
        find_components (B, Map, cn, cnode, Part, Bnz,
            CParent, Cstack, &top, Imap, Common) ;
        /* done using Imap as workspace for Queue ] */
    }
    /* contents of Map [0..cn-1] no longer needed ] */
    }

    /* done using Cmember as workspace for Cmap ] */
    /* done using Perm as workspace for Cstack ] */

    /* ---------------------------------------------------------------------- */
    /* place nodes removed via compression into their proper component */
    /* ---------------------------------------------------------------------- */

    /* At this point, all nodes are of Type 1, 2, or 3, as defined above. */

    for (i = 0 ; i < n ; i++)
    {
    /* find the repnode cnode that contains node i */
    j = FLIP (Flag [i]) ;
    PRINT2 (("\nfind component for "ID", in: "ID"\n", i, j)) ;
    ASSERT (j >= 0 && j < n) ;
    DEBUG (cnt = 0) ;
    while (CParent [j] == -2)
    {
        j = FLIP (Flag [j]) ;
        PRINT2 (("    walk up to "ID" ", j)) ;
        ASSERT (j >= 0 && j < n) ;
        PRINT2 ((" CParent "ID"\n", CParent [j])) ;
        ASSERT (cnt < n) ;
        DEBUG (cnt++) ;
    }
    cnode = j ;
    ASSERT (cnode >= 0 && cnode < n) ;
    ASSERT (CParent [cnode] >= EMPTY && CParent [cnode] < n) ;
    PRINT2 (("i "ID" is in component with cnode "ID"\n", i, cnode)) ;
    ASSERT (Flag [cnode] == FLIP (cnode)) ;

    /* Mark all nodes along the path from i to cnode as being in the
     * component whos repnode is cnode.  Perform path compression.  */
    j = FLIP (Flag [i]) ;
    Flag [i] = FLIP (cnode) ;
    DEBUG (cnt = 0) ;
    while (CParent [j] == -2)
    {
        ASSERT (j >= 0 && j < n) ;
        jnext = FLIP (Flag [j]) ;
        PRINT2 (("    "ID" walk "ID" set cnode to "ID"\n", i, j, cnode)) ;
        ASSERT (cnt < n) ;
        DEBUG (cnt++) ;
        Flag [j] = FLIP (cnode) ;
        j = jnext ;
    }
    }

    /* At this point, all nodes fall into Types 1 or 2, as defined above. */

#ifndef NDEBUG
    for (j = 0 ; j < n ; j++)
    {
    PRINT2 (("j %d CParent %d  ", j, CParent [j])) ;
    if (CParent [j] >= EMPTY && CParent [j] < n)
    {
        /* case 1: j is a repnode of a component */
        cnode = j ;
        PRINT2 ((" a repnode\n")) ;
    }
    else
    {
        /* case 2: j is not a repnode of a component */
        cnode = FLIP (Flag [j]) ;
        PRINT2 ((" repnode is %d\n", cnode)) ;
        ASSERT (cnode >= 0 && cnode < n) ;
        ASSERT (CParent [cnode] >= EMPTY && CParent [cnode] < n) ;
    }
    ASSERT (Flag [cnode] == FLIP (cnode)) ;
    /* case 3 no longer holds */
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* free workspace */
    /* ---------------------------------------------------------------------- */

    C->ncol = n ;   /* restore size for memory usage statistics */
    CHOLMOD(free_sparse) (&C, Common) ;
    CHOLMOD(free_sparse) (&B, Common) ;
    CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
    CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;

    /* ---------------------------------------------------------------------- */
    /* handle dense nodes */
    /* ---------------------------------------------------------------------- */

    /* The separator tree has nodes with either no children or two or more
     * children - with one exception.  There may exist a single root node with
     * exactly one child, which holds the dense rows/columns of the matrix.
     * Delete this node if it exists. */

    if (ndense > 0)
    {
    ASSERT (CParent [cdense] == EMPTY) ;    /* cdense has no parent */
    /* find the children of cdense */
    nchild = 0 ;
    for (j = 0 ; j < n ; j++)
    {
        if (CParent [j] == cdense)
        {
        nchild++ ;
        child = j ;
        }
    }
    if (nchild == 1)
    {
        /* the cdense node has just one child; merge the two nodes */
        PRINT1 (("root has one child\n")) ;
        CParent [cdense] = -2 ;     /* cdense is deleted */
        CParent [child] = EMPTY ;       /* child becomes a root */
        for (j = 0 ; j < n ; j++)
        {
        if (Flag [j] == FLIP (cdense))
        {
            /* j is a dense node */
            PRINT1 (("dense %d\n", j)) ;
            Flag [j] = FLIP (child) ;
        }
        }
    }
    }

    /* ---------------------------------------------------------------------- */
    /* postorder the components */
    /* ---------------------------------------------------------------------- */

    DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) if (CParent [j] != -2) cnt++) ;

    /* use Cmember as workspace for Post [ */
    Post = Cmember ;

    /* cholmod_postorder uses Head and Iwork [0..2n].  It does not use Flag,
     * which here holds the mapping of nodes to repnodes.  It ignores all nodes
     * for which CParent [j] < -1, so it operates just on the repnodes. */
    /* workspace: Head (n), Iwork (2*n) */
    ncomponents = CHOLMOD(postorder) (CParent, n, NULL, Post, Common) ;
    ASSERT (cnt == ncomponents) ;

    /* use Iwork [0..n-1] as workspace for Ipost ( */
    Ipost = Iwork ;
    DEBUG (for (j = 0 ; j < n ; j++) Ipost [j] = EMPTY) ;

    /* compute inverse postorder */
    for (c = 0 ; c < ncomponents ; c++)
    {
    cnode = Post [c] ;
    ASSERT (cnode >= 0 && cnode < n) ;
    Ipost [cnode] = c ;
    ASSERT (Head [c] == EMPTY) ;
    }

    /* adjust the parent array */
    /* Iwork [n..2n-1] used for NewParent [ */
    NewParent = Iwork + n ;
    for (c = 0 ; c < ncomponents ; c++)
    {
    parent = CParent [Post [c]] ;
    NewParent [c] = (parent == EMPTY) ? EMPTY : (Ipost [parent]) ;
    }
    for (c = 0 ; c < ncomponents ; c++)
    {
    CParent [c] = NewParent [c] ;
    }
    ASSERT (CHOLMOD(dump_parent) (CParent, ncomponents, "CParent", Common)) ;

    /* Iwork [n..2n-1] no longer needed for NewParent ] */
    /* Cmember no longer needed for Post ] */

#ifndef NDEBUG
    /* count the number of children of each node */
    for (c = 0 ; c < ncomponents ; c++)
    {
    Cmember [c] = 0 ;
    }
    for (c = 0 ; c < ncomponents ; c++)
    {
    if (CParent [c] != EMPTY) Cmember [CParent [c]]++ ;
    }
    for (c = 0 ; c < ncomponents ; c++)
    {
    /* a node is either a leaf, or has 2 or more children */
    ASSERT (Cmember [c] == 0 || Cmember [c] >= 2) ;
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* place each node in its component */
    /* ---------------------------------------------------------------------- */

    for (j = 0 ; j < n ; j++)
    {
    /* node j is in the cth component, whose repnode is cnode */
    cnode = FLIP (Flag [j]) ;
    PRINT2 (("j "ID"  flag "ID" cnode "ID"\n",
            j, Flag [j], FLIP (Flag [j]))) ;
    ASSERT (cnode >= 0 && cnode < n) ;
    c = Ipost [cnode] ;
    ASSERT (c >= 0 && c < ncomponents) ;
    Cmember [j] = c ;
    }

    /* Flag no longer needed for the node-to-component mapping */

    /* done using Iwork [0..n-1] as workspace for Ipost ) */

    /* ---------------------------------------------------------------------- */
    /* clear the Flag array */
    /* ---------------------------------------------------------------------- */

    Common->mark = EMPTY ;
    CHOLMOD_CLEAR_FLAG (Common) ;
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

    /* ---------------------------------------------------------------------- */
    /* find the permutation */
    /* ---------------------------------------------------------------------- */

    PRINT1 (("nd_camd: %d A->stype %d\n", nd_camd, A->stype)) ;

    if (nd_camd)
    {

    /* ------------------------------------------------------------------ */
    /* apply camd, csymamd, or ccolamd using the Cmember constraints */
    /* ------------------------------------------------------------------ */

    if (A->stype != 0)
    {
        /* ordering A+A', so fset and fsize are ignored.
         * Add the upper/lower part to a symmetric lower/upper matrix by
         * converting to unsymmetric mode
         * workspace: Iwork (nrow) */
        B = CHOLMOD(copy) (A, 0, -1, Common) ;
        if (Common->status < CHOLMOD_OK)
        {
        PRINT0 (("make symmetric failed\n")) ;
        return (EMPTY) ;
        }
        ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;
        PRINT2 (("nested dissection (2)\n")) ;
        B->stype = -1 ;
        if (nd_camd == 2)
        {
        /* workspace:  Head (nrow+1), Iwork (nrow) if symmetric-upper */
        ok = CHOLMOD(csymamd) (B, Cmember, Perm, Common) ;
        }
        else
        {
        /* workspace: Head (nrow), Iwork (4*nrow) */
        ok = CHOLMOD(camd) (B, NULL, 0, Cmember, Perm, Common) ;
        }
        CHOLMOD(free_sparse) (&B, Common) ;
        if (!ok)
        {
        /* failed */
        PRINT0 (("camd/csymamd failed\n")) ;
        return (EMPTY) ;
        }
    }
    else
    {
        /* ordering A*A' or A(:,f)*A(:,f)' */
        /* workspace: Iwork (nrow if no fset; MAX(nrow,ncol) if fset) */
        if (!CHOLMOD(ccolamd) (A, fset, fsize, Cmember, Perm, Common))
        {
        /* ccolamd failed */
        PRINT2 (("ccolamd failed\n")) ;
        return (EMPTY) ;
        }
    }

    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* natural ordering of each component */
    /* ------------------------------------------------------------------ */

    /* use Iwork [0..n-1] for Next [ */
    Next = Iwork  ;

    /* ------------------------------------------------------------------ */
    /* place the nodes in link lists, one list per component */
    /* ------------------------------------------------------------------ */

    /* do so in reverse order, to preserve original ordering */
    for (j = n-1 ; j >= 0 ; j--)
    {
        /* node j is in the cth component */
        c = Cmember [j] ;
        ASSERT (c >= 0 && c < ncomponents) ;
        /* place node j in link list for component c */
        Next [j] = Head [c] ;
        Head [c] = j ;
    }

    /* ------------------------------------------------------------------ */
    /* order each node in each component */
    /* ------------------------------------------------------------------ */

    k = 0 ;
    for (c = 0 ; c < ncomponents ; c++)
    {
        for (j = Head [c] ; j != EMPTY ; j = Next [j])
        {
        Perm [k++] = j ;
        }
        Head [c] = EMPTY ;
    }
    ASSERT (k == n) ;

    /* done using Iwork [0..n-1] for Next ] */
    }

    /* ---------------------------------------------------------------------- */
    /* clear workspace and return number of components */
    /* ---------------------------------------------------------------------- */

    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
    return (ncomponents) ;
}

/* ========================================================================== */
/* === cholmod_collapse_septree ============================================= */
/* ========================================================================== */

/* cholmod_nested_dissection returns the separator tree that was used in the
 * constrained minimum degree algorithm.  Parameter settings (nd_small,
 * nd_oksep, etc) that give a good fill-reducing ordering may give too fine of
 * a separator tree for other uses (parallelism, multi-level LPDASA, etc).  This
 * function takes as input the separator tree computed by
 * cholmod_nested_dissection, and collapses selected subtrees into single
 * nodes.  A subtree is collapsed if its root node (the separator) is large
 * compared to the total number of nodes in the subtree, or if the subtree is
 * small.  Note that the separator tree may actually be a forest.
 *
 * nd_oksep and nd_small act just like the ordering parameters in Common.
 * Returns the new number of nodes in the separator tree.
 */

UF_long CHOLMOD(collapse_septree)
(
    /* ---- input ---- */
    size_t n,       /* # of nodes in the graph */
    size_t ncomponents, /* # of nodes in the separator tree (must be <= n) */
    double nd_oksep,    /* collapse if #sep >= nd_oksep * #nodes in subtree */
    size_t nd_small,    /* collapse if #nodes in subtree < nd_small */
    /* ---- in/out --- */
    Int *CParent,   /* size ncomponents; from cholmod_nested_dissection */
    Int *Cmember,   /* size n; from cholmod_nested_dissection */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *First, *Count, *Csubtree, *W, *Map ;
    Int c, j, k, nc, sepsize, total_weight, parent, nc_new, first ;
    int collapse = FALSE, ok = TRUE ;
    size_t s ;

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (EMPTY) ;
    RETURN_IF_NULL (CParent, EMPTY) ;
    RETURN_IF_NULL (Cmember, EMPTY) ;
    if (n < ncomponents)
    {
    ERROR (CHOLMOD_INVALID, "invalid separator tree") ;
    return (EMPTY) ;
    }
    Common->status = CHOLMOD_OK ;
    nc = ncomponents ;
    if (n <= 1 || ncomponents <= 1)
    {
    /* no change; tree is one node already */
    return (nc) ;
    }

    nd_oksep = MAX (0, nd_oksep) ;
    nd_oksep = MIN (1, nd_oksep) ;
    nd_small = MAX (4, nd_small) ;

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    /* s = 3*ncomponents */
    s = CHOLMOD(mult_size_t) (ncomponents, 3, &ok) ;
    if (!ok)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    return (EMPTY) ;
    }
    CHOLMOD(allocate_work) (0, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (EMPTY) ;
    }
    W = Common->Iwork ;
    Count    = W ; W += ncomponents ;       /* size ncomponents */
    Csubtree = W ; W += ncomponents ;       /* size ncomponents */
    First    = W ; W += ncomponents ;       /* size ncomponents */

    /* ---------------------------------------------------------------------- */
    /* find the first descendant of each node of the separator tree */
    /* ---------------------------------------------------------------------- */

    for (c = 0 ; c < nc ; c++)
    {
    First [c] = EMPTY ;
    }
    for (k = 0 ; k < nc ; k++)
    {
    for (c = k ; c != EMPTY && First [c] == -1 ; c = CParent [c])
    {
        ASSERT (c >= 0 && c < nc) ;
        First [c] = k ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* find the number of nodes of the graph in each node of the tree */
    /* ---------------------------------------------------------------------- */

    for (c = 0 ; c < nc ; c++)
    {
    Count [c] = 0 ;
    }
    for (j = 0 ; j < (Int) n ; j++)
    {
    ASSERT (Cmember [j] >= 0 && Cmember [j] < nc) ;
    Count [Cmember [j]]++ ;
    }

    /* ---------------------------------------------------------------------- */
    /* find the number of nodes in each subtree */
    /* ---------------------------------------------------------------------- */

    for (c = 0 ; c < nc ; c++)
    {
    /* each subtree includes its root */
    Csubtree [c] = Count [c] ;
    PRINT1 ((ID" size "ID" parent "ID" first "ID"\n",
        c, Count [c], CParent [c], First [c])) ;
    }

    for (c = 0 ; c < nc ; c++)
    {
    /* add the subtree of the child, c, into the count of its parent */
    parent = CParent [c] ;
    ASSERT (parent >= EMPTY && parent < nc) ;
    if (parent != EMPTY)
    {
        Csubtree [parent] += Csubtree [c] ;
    }
    }

#ifndef NDEBUG
    /* the sum of the roots should be n */
    j = 0 ;
    for (c = 0 ; c < nc ; c++) if (CParent [c] == EMPTY) j += Csubtree [c] ;
    ASSERT (j == (Int) n) ;
#endif

    /* ---------------------------------------------------------------------- */
    /* find subtrees to collapse */
    /* ---------------------------------------------------------------------- */

    /* consider all nodes in reverse post-order */
    for (c = nc-1 ; c >= 0 ; c--)
    {
    /* consider the subtree rooted at node c */
    sepsize = Count [c] ;
    total_weight = Csubtree [c] ;
    PRINT1 (("Node "ID" sepsize "ID" subtree "ID" ratio %g\n", c, sepsize,
        total_weight, ((double) sepsize)/((double) total_weight))) ;
    first = First [c] ;
    if (first < c &&    /* c must not be a leaf */
       (sepsize > nd_oksep * total_weight || total_weight < (int) nd_small))
    {
        /* this separator is too large, or the subtree is too small.
         * collapse the tree, by converting the entire subtree rooted at
         * c into a single node.  The subtree consists of all nodes from
         * First[c] to the root c.  Flag all nodes from First[c] to c-1
         * as dead.
         */
        collapse = TRUE ;
        for (k = first ; k < c ; k++)
        {
        CParent [k] = -2 ;
        PRINT1 (("   collapse node "ID"\n", k)) ;
        }
        /* continue at the next node, first-1 */
        c = first ;
    }
    }

    PRINT1 (("collapse: %d\n", collapse)) ;

    /* ---------------------------------------------------------------------- */
    /* compress the tree */
    /* ---------------------------------------------------------------------- */

    Map = Count ;   /* Count no longer needed */

    nc_new = nc ;
    if (collapse)
    {
    nc_new = 0 ;
    for (c = 0 ; c < nc ; c++)
    {
        Map [c] = nc_new ;
        if (CParent [c] >= EMPTY)
        {
        /* node c is alive, and becomes node Map[c] in the new tree.
         * Increment nc_new for the next node c. */
        nc_new++ ;
        }
    }
    PRINT1 (("Collapse the tree from "ID" to "ID" nodes\n", nc, nc_new)) ;
    ASSERT (nc_new > 0) ;
    for (c = 0 ; c < nc ; c++)
    {
        parent = CParent [c] ;
        if (parent >= EMPTY)
        {
        /* node c is alive */
        CParent [Map [c]] = (parent == EMPTY) ? EMPTY : Map [parent] ;
        }
    }
    for (j = 0 ; j < (Int) n ; j++)
    {
        PRINT1 (("j "ID" Cmember[j] "ID" Map[Cmember[j]] "ID"\n",
        j, Cmember [j], Map [Cmember [j]])) ;
        Cmember [j] = Map [Cmember [j]] ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* return new size of separator tree */
    /* ---------------------------------------------------------------------- */

    return (nc_new) ;
}
#endif