amesos_cholmod_nesdis.c

Go to the documentation of this file.

/* ========================================================================== */
/* === 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 "amesos_cholmod_internal.h"
#include "amesos_cholmod_partition.h"
#include "amesos_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.
 *
 * 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.
 *
 * workspace: Flag (nrow)
 */

static UF_long clear_flag (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 ;
  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 ;
    Int *Bp, *Bi, *Flag ;

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

    PRINT2 (("find components: cn %d\n", cn)) ;
    Flag = Common->Flag ;     /* size n */
    Common->mark = EMPTY ;      /* force initialization of Flag array */
    mark = clear_flag (Common) ;    /* clear Flag but preserve Flag [i]<EMPTY */
    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 */
    sj = 0 ;
    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 ;
      }
  }
    }

    /* clear Flag array, but preserve Flag [i] < EMPTY */
    clear_flag (Common) ;
}


/* ========================================================================== */
/* === 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 */
    /* ---------------------------------------------------------------------- */

    ASSERT (CHOLMOD(dump_sparse) (A, "A for NESDIS:", Common) >= 0) ;

    /* 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 (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