amesos_btf_strongcomp.c

Go to the documentation of this file.

/* ========================================================================== */
/* === BTF_STRONGCOMP ======================================================= */
/* ========================================================================== */

/* Finds the strongly connected components of a graph, or equivalently, permutes
 * the matrix into upper block triangular form.  See btf.h for more details.
 * Input matrix and Q are not checked on input.
 *
 * Copyright (c) 2004-2007.  Tim Davis, University of Florida,
 * with support from Sandia National Laboratories.  All Rights Reserved.
 */

#include "amesos_btf_decl.h"
#include "amesos_btf_internal.h"

#define UNVISITED (-2)      /* Flag [j] = UNVISITED if node j not visited yet */
#define UNASSIGNED (-1)     /* Flag [j] = UNASSIGNED if node j has been visited,
           * but not yet assigned to a strongly-connected
           * component (aka block).  Flag [j] = k (k in the
           * range 0 to nblocks-1) if node j has been visited
           * (and completed, with its postwork done) and
           * assigned to component k. */

/* This file contains two versions of the depth-first-search, a recursive one
 * and a non-recursive one.  By default, the non-recursive one is used. */

#ifndef RECURSIVE

/* ========================================================================== */
/* === dfs: non-recursive version (default) ================================= */
/* ========================================================================== */

/* Perform a depth-first-search of a graph, stored in an adjacency-list form.
 * The row indices of column j (equivalently, the out-adjacency list of node j)
 * are stored in Ai [Ap[j] ... Ap[j+1]-1].  Self-edge (diagonal entries) are
 * ignored.  Ap[0] must be zero, and thus nz = Ap[n] is the number of entries
 * in the matrix (or edges in the graph).  The row indices in each column need
 * not be in any particular order.  If an input column permutation is given,
 * node j (in the permuted matrix A*Q) is located in
 * Ai [Ap[Q[j]] ... Ap[Q[j]+1]-1].  This Q can be the same as the Match array
 * output from the maxtrans routine, for a square matrix that is structurally
 * full rank.
 *
 * The algorithm is from the paper by Robert E. Tarjan, "Depth-first search and
 * linear graph algorithms," SIAM Journal on Computing, vol. 1, no. 2,
 * pp. 146-160, 1972.  The time taken by strongcomp is O(nnz(A)).
 *
 * See also MC13A/B in the Harwell subroutine library (Iain S. Duff and John
 * K. Reid, "Algorithm 529: permutations to block triangular form," ACM Trans.
 * on Mathematical Software, vol. 4, no. 2, pp. 189-192, 1978, and "An
 * implementation of Tarjan's algorithm for the block triangular form of a
 * matrix," same journal, pp. 137-147.  This code is implements the same
 * algorithm as MC13A/B, except that the data structures are very different.
 * Also, unlike MC13A/B, the output permutation preserves the natural ordering
 * within each block.
 */

static void amesos_dfs
(
    /* inputs, not modified on output: */
    Int j,    /* start the DFS at node j */
    Int Ap [ ],   /* size n+1, column pointers for the matrix A */
    Int Ai [ ],   /* row indices, size nz = Ap [n] */
    Int Q [ ],    /* input column permutation */

    /* inputs, modified on output (each array is of size n): */
    Int Time [ ], /* Time [j] = "time" that node j was first visited */
    Int Flag [ ], /* Flag [j]: see above */
    Int Low [ ],  /* Low [j]: see definition below */
    Int *p_nblocks, /* number of blocks (aka strongly-connected-comp.)*/
    Int *p_timestamp, /* current "time" */

    /* workspace, not defined on input or output: */
    Int Cstack [ ], /* size n, output stack to hold nodes of components */
    Int Jstack [ ], /* size n, stack for the variable j */
    Int Pstack [ ]  /* size n, stack for the variable p */
)
{
    /* ---------------------------------------------------------------------- */
    /* local variables, and initializations */
    /* ---------------------------------------------------------------------- */

    /* local variables, but "global" to all DFS levels: */
    Int chead ;     /* top of Cstack */
    Int jhead ;     /* top of Jstack and Pstack */

    /* variables that are purely local to any one DFS level: */
    Int i ;     /* edge (j,i) considered; i can be next node to traverse */
    Int parent ;    /* parent of node j in the DFS tree */
    Int pend ;      /* one past the end of the adjacency list for node j */
    Int jj ;      /* column j of A*Q is column jj of the input matrix A */

    /* variables that need to be pushed then popped from the stack: */
    Int p ;     /* current index into the adj. list for node j */
    /* the variables j and p are stacked in Jstack and Pstack */

    /* local copies of variables in the calling routine */
    Int nblocks   = *p_nblocks ;
    Int timestamp = *p_timestamp ;

    /* ---------------------------------------------------------------------- */
    /* start a DFS at node j (same as the recursive call dfs (EMPTY, j)) */
    /* ---------------------------------------------------------------------- */

    chead = 0 ;       /* component stack is empty */
    jhead = 0 ;       /* Jstack and Pstack are empty */
    Jstack [0] = j ;      /* put the first node j on the Jstack */
    ASSERT (Flag [j] == UNVISITED) ;

    while (jhead >= 0)
    {
  j = Jstack [jhead] ;      /* grab the node j from the top of Jstack */

  /* determine which column jj of the A is column j of A*Q */
  jj = (Q == (Int *) NULL) ? (j) : (BTF_UNFLIP (Q [j])) ;
  pend = Ap [jj+1] ;      /* j's row index list ends at Ai [pend-1] */

  if (Flag [j] == UNVISITED)
  {

      /* -------------------------------------------------------------- */
      /* prework at node j */
      /* -------------------------------------------------------------- */

      /* node j is being visited for the first time */
      Cstack [++chead] = j ;      /* push j onto the stack */
      timestamp++ ;       /* get a timestamp */
      Time [j] = timestamp ;      /* give the timestamp to node j */
      Low [j] = timestamp ;
      Flag [j] = UNASSIGNED ;     /* flag node j as visited */

      /* -------------------------------------------------------------- */
      /* set Pstack [jhead] to the first entry in column j to scan */
      /* -------------------------------------------------------------- */

      Pstack [jhead] = Ap [jj] ;
  }

  /* ------------------------------------------------------------------ */
  /* DFS rooted at node j (start it, or continue where left off) */
  /* ------------------------------------------------------------------ */

  for (p = Pstack [jhead] ; p < pend ; p++)
  {
      i = Ai [p] ;    /* examine the edge from node j to node i */
      if (Flag [i] == UNVISITED)
      {
    /* Node i has not been visited - start a DFS at node i.
     * Keep track of where we left off in the scan of adjacency list
     * of node j so we can restart j where we left off. */
    Pstack [jhead] = p + 1 ;
    /* Push i onto the stack and immediately break
     * so we can recurse on node i. */
    Jstack [++jhead] = i ;
    ASSERT (Time [i] == EMPTY) ;
    ASSERT (Low [i] == EMPTY) ;
    /* break here to do what the recursive call dfs (j,i) does */
    break ;
      }
      else if (Flag [i] == UNASSIGNED)
      {
    /* Node i has been visited, but still unassigned to a block
     * this is a back or cross edge if Time [i] < Time [j].
     * Note that i might equal j, in which case this code does
     * nothing. */
    ASSERT (Time [i] > 0) ;
    ASSERT (Low [i] > 0) ;
    Low [j] = MIN (Low [j], Time [i]) ;
      }
  }

  if (p == pend)
  {
      /* If all adjacent nodes of j are already visited, pop j from
       * Jstack and do the post work for node j.  This also pops p
       * from the Pstack. */
      jhead-- ;

      /* -------------------------------------------------------------- */
      /* postwork at node j */
      /* -------------------------------------------------------------- */

      /* determine if node j is the head of a component */
      if (Low [j] == Time [j])
      {
    /* pop all nodes in this SCC from Cstack */
    while (TRUE)
    {
        ASSERT (chead >= 0) ; /* stack not empty (j in it) */
        i = Cstack [chead--] ;  /* pop a node from the Cstack */
        ASSERT (i >= 0) ;
        ASSERT (Flag [i] == UNASSIGNED) ;
        Flag [i] = nblocks ;  /* assign i to current block */
        if (i == j) break ;   /* current block ends at j */
    }
    nblocks++ ; /* one more block has been found */
      }
      /* update Low [parent], if the parent exists */
      if (jhead >= 0)
      {
    parent = Jstack [jhead] ;
    Low [parent] = MIN (Low [parent], Low [j]) ;
      }
  }
    }

    /* ---------------------------------------------------------------------- */
    /* cleanup: update timestamp and nblocks */
    /* ---------------------------------------------------------------------- */

    *p_timestamp = timestamp ;
    *p_nblocks   = nblocks ;
}

#else

/* ========================================================================== */
/* === dfs: recursive version (only for illustration) ======================= */
/* ========================================================================== */

/* The following is a recursive version of dfs, which computes identical results
 * as the non-recursive dfs.  It is included here because it is easier to read.
 * Compare the comments in the code below with the identical comments in the
 * non-recursive code above, and that will help you see the correlation between
 * the two routines.
 *
 * This routine can cause stack overflow, and is thus not recommended for heavy
 * usage, particularly for large matrices.  To help in delaying stack overflow,
 * global variables are used, reducing the amount of information each call to
 * dfs places on the call/return stack (the integers i, j, p, parent, and the
 * return address).  Note that this means the recursive code is not thread-safe.
 * To try this version, compile the code with -DRECURSIVE or include the
 * following line at the top of this file:

#define RECURSIVE

 */

static Int  /* for recursive illustration only, not for production use */
    chead, timestamp, nblocks, n, *Ap, *Ai, *Flag, *Cstack, *Time, *Low,
    *P, *R, *Q ;

static void amesos_dfs
(
    Int parent,   /* came from parent node */
    Int j   /* at node j in the DFS */
)
{
    Int p ;     /* current index into the adj. list for node j */
    Int i ;     /* edge (j,i) considered; i can be next node to traverse */
    Int jj ;      /* column j of A*Q is column jj of the input matrix A */

    /* ---------------------------------------------------------------------- */
    /* prework at node j */
    /* ---------------------------------------------------------------------- */

    /* node j is being visited for the first time */
    Cstack [++chead] = j ;      /* push j onto the stack */
    timestamp++ ;       /* get a timestamp */
    Time [j] = timestamp ;      /* give the timestamp to node j */
    Low [j] = timestamp ;
    Flag [j] = UNASSIGNED ;     /* flag node j as visited */

    /* ---------------------------------------------------------------------- */
    /* DFS rooted at node j */
    /* ---------------------------------------------------------------------- */

    /* determine which column jj of the A is column j of A*Q */
    jj = (Q == (Int *) NULL) ? (j) : (BTF_UNFLIP (Q [j])) ;
    for (p = Ap [jj] ; p < Ap [jj+1] ; p++)
    {
  i = Ai [p] ;    /* examine the edge from node j to node i */
  if (Flag [i] == UNVISITED)
  {
      /* Node i has not been visited - start a DFS at node i. */
      amesos_dfs (j, i) ;
  }
  else if (Flag [i] == UNASSIGNED)
  {
      /* Node i has been visited, but still unassigned to a block
       * this is a back or cross edge if Time [i] < Time [j].
       * Note that i might equal j, in which case this code does
       * nothing. */
      Low [j] = MIN (Low [j], Time [i]) ;
  }
    }

    /* ---------------------------------------------------------------------- */
    /* postwork at node j */
    /* ---------------------------------------------------------------------- */

    /* determine if node j is the head of a component */
    if (Low [j] == Time [j])
    {
  /* pop all nodes in this strongly connected component from Cstack */
  while (TRUE)
  {
      i = Cstack [chead--] ;  /* pop a node from the Cstack */
      Flag [i] = nblocks ;  /* assign node i to current block */
      if (i == j) break ;   /* current block ends at node j */
  }
  nblocks++ ; /* one more block has been found */
    }
    /* update Low [parent] */
    if (parent != EMPTY)
    {
  /* Note that this could be done with Low[j] = MIN(Low[j],Low[i]) just
   * after the dfs (j,i) statement above, and then parent would not have
   * to be an input argument.  Putting it here places all the postwork
   * for node j in one place, thus making the non-recursive DFS easier. */
  Low [parent] = MIN (Low [parent], Low [j]) ;
    }
}

#endif

/* ========================================================================== */
/* === btf_strongcomp ======================================================= */
/* ========================================================================== */

#ifndef RECURSIVE

Int BTF(strongcomp) /* return # of strongly connected components */
(
    /* input, not modified: */
    Int n,      /* A is n-by-n in compressed column form */
    Int Ap [ ],     /* size n+1 */
    Int Ai [ ],     /* size nz = Ap [n] */

    /* optional input, modified (if present) on output: */
    Int Q [ ],      /* size n, input column permutation.  The permutation Q can
         * include a flag which indicates an unmatched row.
         * jold = BTF_UNFLIP (Q [jnew]) is the permutation;
         * this function ingnores these flags.  On output, it is
         * modified according to the permutation P. */

    /* output, not defined on input: */
    Int P [ ],      /* size n.  P [k] = j if row and column j are kth row/col
         * in permuted matrix. */
    Int R [ ],      /* size n+1.  kth block is in rows/cols R[k] ... R[k+1]-1
         * of the permuted matrix. */

    /* workspace, not defined on input or output: */
    Int Work [ ]    /* size 4n */
)

#else

Int BTF(strongcomp) /* recursive version - same as above except for Work size */
(
    Int n_in,
    Int Ap_in [ ],
    Int Ai_in [ ],
    Int Q_in [ ],
    Int P_in [ ],
    Int R_in [ ],
    Int Work [ ]    /* size 2n */
)

#endif

{
    Int j, k, b ;

#ifndef RECURSIVE
    Int timestamp, nblocks, *Flag, *Cstack, *Time, *Low, *Jstack, *Pstack ;
#else
    n = n_in ;
    Ap = Ap_in ;
    Ai = Ai_in ;
    Q = Q_in ;
    P = P_in ;
    R = R_in ;
    chead = EMPTY ;
#endif

    /* ---------------------------------------------------------------------- */
    /* get and initialize workspace */
    /* ---------------------------------------------------------------------- */

    /* timestamp is incremented each time a new node is visited.
     *
     * Time [j] is the timestamp given to node j.
     *
     * Low [j] is the lowest timestamp of any node reachable from j via either
     * a path to any descendent of j in the DFS tree, or via a single edge to
     * an either an ancestor (a back edge) or another node that's neither an
     * ancestor nor a descendant (a cross edge).  If Low [j] is equal to
     * the timestamp of node j (Time [j]), then node j is the "head" of a
     * strongly connected component (SCC).  That is, it is the first node
     * visited in its strongly connected component, and the DFS subtree rooted
     * at node j spans all the nodes of the strongly connected component.
     *
     * The term "block" and "component" are used interchangebly in this code;
     * "block" being a matrix term and "component" being a graph term for the
     * same thing.
     *
     * When a node is visited, it is placed on the Cstack (for "component"
     * stack).  When node j is found to be an SCC head, all the nodes from the
     * top of the stack to node j itself form the nodes in the SCC.  This Cstack
     * is used for both the recursive and non-recursive versions.
     */

    Time   = Work ; Work += n ;
    Flag   = Work ; Work += n ;
    Low    = P ;    /* use output array P as workspace for Low */
    Cstack = R ;    /* use output array R as workspace for Cstack */

#ifndef RECURSIVE
    /* stack for non-recursive dfs */
    Jstack = Work ; Work += n ;     /* stack for j */
    Pstack = Work ;       /* stack for p */
#endif

    for (j = 0 ; j < n ; j++)
    {
  Flag [j] = UNVISITED ;
  Low [j] = EMPTY ;
  Time [j] = EMPTY ;
#ifndef NDEBUG
  Cstack [j] = EMPTY ;
#ifndef RECURSIVE
  Jstack [j] = EMPTY ;
  Pstack [j] = EMPTY ;
#endif
#endif
    }

    timestamp = 0 ; /* each node given a timestamp when it is visited */
    nblocks = 0 ; /* number of blocks found so far */

    /* ---------------------------------------------------------------------- */
    /* find the connected components via a depth-first-search */
    /* ---------------------------------------------------------------------- */

    for (j = 0 ; j < n ; j++)
    {
  /* node j is unvisited or assigned to a block. Cstack is empty. */
  ASSERT (Flag [j] == UNVISITED || (Flag [j] >= 0 && Flag [j] < nblocks));
  if (Flag [j] == UNVISITED)
  {
#ifndef RECURSIVE
      /* non-recursive dfs (default) */
      amesos_dfs (j, Ap, Ai, Q, Time, Flag, Low, &nblocks, &timestamp,
        Cstack, Jstack, Pstack) ;
#else
      /* recursive dfs (for illustration only) */
      ASSERT (chead == EMPTY) ;
      amesos_dfs (EMPTY, j) ;
      ASSERT (chead == EMPTY) ;
#endif
  }
    }
    ASSERT (timestamp == n) ;

    /* ---------------------------------------------------------------------- */
    /* construct the block boundary array, R */
    /* ---------------------------------------------------------------------- */

    for (b = 0 ; b < nblocks ; b++)
    {
  R [b] = 0 ;
    }
    for (j = 0 ; j < n ; j++)
    {
  /* node j has been assigned to block b = Flag [j] */
  ASSERT (Time [j] > 0 && Time [j] <= n) ;
  ASSERT (Low [j] > 0 && Low [j] <= n) ;
  ASSERT (Flag [j] >= 0 && Flag [j] < nblocks) ;
  R [Flag [j]]++ ;
    }
    /* R [b] is now the number of nodes in block b.  Compute cumulative sum
     * of R, using Time [0 ... nblocks-1] as workspace. */
    Time [0] = 0 ;
    for (b = 1 ; b < nblocks ; b++)
    {
  Time [b] = Time [b-1] + R [b-1] ;
    }
    for (b = 0 ; b < nblocks ; b++)
    {
  R [b] = Time [b] ;
    }
    R [nblocks] = n ;

    /* ---------------------------------------------------------------------- */
    /* construct the permutation, preserving the natural order */
    /* ---------------------------------------------------------------------- */

#ifndef NDEBUG
    for (k = 0 ; k < n ; k++)
    {
  P [k] = EMPTY ;
    }
#endif

    for (j = 0 ; j < n ; j++)
    {
  /* place column j in the permutation */
  P [Time [Flag [j]]++] = j ;
    }

#ifndef NDEBUG
    for (k = 0 ; k < n ; k++)
    {
  ASSERT (P [k] != EMPTY) ;
    }
#endif

    /* Now block b consists of the nodes k1 to k2-1 in the permuted matrix,
     * where k1 = R [b] and k2 = R [b+1].  Row and column j of the original
     * matrix becomes row and column P [k] of the permuted matrix.  The set of
     * of rows/columns (nodes) in block b is given by P [k1 ... k2-1], and this
     * set is sorted in ascending order.  Thus, if the matrix consists of just
     * one block, P is the identity permutation. */

    /* ---------------------------------------------------------------------- */
    /* if Q is present on input, set Q = Q*P' */
    /* ---------------------------------------------------------------------- */

    if (Q != (Int *) NULL)
    {
  /* We found a symmetric permutation P for the matrix A*Q.  The overall
   * permutation is thus P*(A*Q)*P'.  Set Q=Q*P' so that the final
   * permutation is P*A*Q.  Use Time as workspace.  Note that this
   * preserves the negative values of Q if the matrix is structurally
   * singular. */
  for (k = 0 ; k < n ; k++)
  {
      Time [k] = Q [P [k]] ;
  }
  for (k = 0 ; k < n ; k++)
  {
      Q [k] = Time [k] ;
  }
    }

    /* ---------------------------------------------------------------------- */
    /* how to traverse the permuted matrix */
    /* ---------------------------------------------------------------------- */

    /* If Q is not present, the following code can be used to traverse the
     * permuted matrix P*A*P'
     *
     *      // compute the inverse of P
     *      for (knew = 0 ; knew < n ; knew++)
     *      {
     *    // row and column kold in the old matrix is row/column knew
     *    // in the permuted matrix P*A*P'
     *    kold = P [knew] ;
     *    Pinv [kold] = knew ;
     *      }
     *      for (b = 0 ; b < nblocks ; b++)
     *      {
     *    // traverse block b of the permuted matrix P*A*P'
     *    k1 = R [b] ;
     *    k2 = R [b+1] ;
     *    nk = k2 - k1 ;
     *    for (jnew = k1 ; jnew < k2 ; jnew++)
     *    {
     *        jold = P [jnew] ;
     *        for (p = Ap [jold] ; p < Ap [jold+1] ; p++)
     *        {
     *      iold = Ai [p] ;
     *      inew = Pinv [iold] ;
     *      // Entry in the old matrix is A (iold, jold), and its
     *      // position in the new matrix P*A*P' is (inew, jnew).
     *      // Let B be the bth diagonal block of the permuted
     *      // matrix.  If inew >= k1, then this entry is in row/
     *      // column (inew-k1, jnew-k1) of the nk-by-nk matrix B.
     *      // Otherwise, the entry is in the upper block triangular
     *      // part, not in any diagonal block.
     *        }
     *    }
     *      }
     *
     * If Q is present replace the above statement
     *    jold = P [jnew] ;
     * with
     *    jold = Q [jnew] ;
     * or
     *    jold = BTF_UNFLIP (Q [jnew]) ;
     *
     * then entry A (iold,jold) in the old (unpermuted) matrix is at (inew,jnew)
     * in the permuted matrix P*A*Q.  Everything else remains the same as the
     * above (simply replace P*A*P' with P*A*Q in the above comments).
     */

    /* ---------------------------------------------------------------------- */
    /* return # of blocks / # of strongly connected components */
    /* ---------------------------------------------------------------------- */

    return (nblocks) ;
}