text: #include "cxs.h"
  /* column counts of LL'=A or LL'=A'A, given parent & post ordering */
  #define HEAD(k,j) (ata ? head [k] : j)
  #define NEXT(J)   (ata ? next [J] : -1)
  static void init_ata (cxs *AT, const CXS_INT *post, CXS_INT *w, CXS_INT **head, CXS_INT **next)
  {
      CXS_INT i, k, p, m = AT->n, n = AT->m, *ATp = AT->p, *ATi = AT->i ;
      *head = w+4*n, *next = w+5*n+1 ;
      for (k = 0 ; k < n ; k++) w [post [k]] = k ;    /* invert post */
      for (i = 0 ; i < m ; i++)
      {
          for (k = n, p = ATp[i] ; p < ATp[i+1] ; p++) k = CXS_MIN (k, w [ATi[p]]);
          (*next) [i] = (*head) [k] ;     /* place row i in linked list k */
          (*head) [k] = i ;
      }
  }
  CXS_INT *cxs_counts (const cxs *A, const CXS_INT *parent, const CXS_INT *post, CXS_INT ata)
  {
      CXS_INT i, j, k, n, m, J, s, p, q, jleaf, *ATp, *ATi, *maxfirst, *prevleaf,
          *ancestor, *head = NULL, *next = NULL, *colcount, *w, *first, *delta ;
      cxs *AT ;
      if (!CXS_CSC (A) || !parent || !post) return (NULL) ;    /* check inputs */
      m = A->m ; n = A->n ;
      s = 4*n + (ata ? (n+m+1) : 0) ;
      delta = colcount = cxs_malloc (n, sizeof (CXS_INT)) ;    /* allocate result */
      w = cxs_malloc (s, sizeof (CXS_INT)) ;                   /* get workspace */
      AT = cxs_transpose (A, 0) ;                          /* AT = A' */
      if (!AT || !colcount || !w) return (cxs_idone (colcount, AT, w, 0)) ;
      ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ; first = w+3*n ;
      for (k = 0 ; k < s ; k++) w [k] = -1 ;      /* clear workspace w [0..s-1] */
      for (k = 0 ; k < n ; k++)                   /* find first [j] */
      {
          j = post [k] ;
          delta [j] = (first [j] == -1) ? 1 : 0 ;  /* delta[j]=1 if j is a leaf */
          for ( ; j != -1 && first [j] == -1 ; j = parent [j]) first [j] = k ;
      }
      ATp = AT->p ; ATi = AT->i ;
      if (ata) init_ata (AT, post, w, &head, &next) ;
      for (i = 0 ; i < n ; i++) ancestor [i] = i ; /* each node in its own set */
      for (k = 0 ; k < n ; k++)
      {
          j = post [k] ;          /* j is the kth node in postordered etree */
          if (parent [j] != -1) delta [parent [j]]-- ;    /* j is not a root */
          for (J = HEAD (k,j) ; J != -1 ; J = NEXT (J))   /* J=j for LL'=A case */
          {
              for (p = ATp [J] ; p < ATp [J+1] ; p++)
              {
                  i = ATi [p] ;
                  q = cxs_leaf (i, j, first, maxfirst, prevleaf, ancestor, &jleaf);
                  if (jleaf >= 1) delta [j]++ ;   /* A(i,j) is in skeleton */
                  if (jleaf == 2) delta [q]-- ;   /* account for overlap in q */
              }
          }
          if (parent [j] != -1) ancestor [j] = parent [j] ;
      }
      for (j = 0 ; j < n ; j++)           /* sum up delta's of each child */
      {
          if (parent [j] != -1) colcount [parent [j]] += colcount [j] ;
      }
      return (cxs_idone (colcount, AT, w, 1)) ;    /* success: free workspace */
  }
  