Amesos Package Browser (Single Doxygen Collection) Development
amesos_btf_l_strongcomp.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === BTF_STRONGCOMP ======================================================= */
00003 /* ========================================================================== */
00004 
00005 /* Finds the strongly connected components of a graph, or equivalently, permutes
00006  * the matrix into upper block triangular form.  See btf.h for more details.
00007  * Input matrix and Q are not checked on input.
00008  *
00009  * Copyright (c) 2004-2007.  Tim Davis, University of Florida,
00010  * with support from Sandia National Laboratories.  All Rights Reserved.
00011  */
00012 
00013 /* This file should make the long int version of BTF */
00014 #define DLONG 1
00015 
00016 #include "amesos_btf_decl.h"
00017 #include "amesos_btf_internal.h"
00018 
00019 #define UNVISITED (-2)      /* Flag [j] = UNVISITED if node j not visited yet */
00020 #define UNASSIGNED (-1)     /* Flag [j] = UNASSIGNED if node j has been visited,
00021            * but not yet assigned to a strongly-connected
00022            * component (aka block).  Flag [j] = k (k in the
00023            * range 0 to nblocks-1) if node j has been visited
00024            * (and completed, with its postwork done) and
00025            * assigned to component k. */
00026 
00027 /* This file contains two versions of the depth-first-search, a recursive one
00028  * and a non-recursive one.  By default, the non-recursive one is used. */
00029 
00030 #ifndef RECURSIVE
00031 
00032 /* ========================================================================== */
00033 /* === dfs: non-recursive version (default) ================================= */
00034 /* ========================================================================== */
00035 
00036 /* Perform a depth-first-search of a graph, stored in an adjacency-list form.
00037  * The row indices of column j (equivalently, the out-adjacency list of node j)
00038  * are stored in Ai [Ap[j] ... Ap[j+1]-1].  Self-edge (diagonal entries) are
00039  * ignored.  Ap[0] must be zero, and thus nz = Ap[n] is the number of entries
00040  * in the matrix (or edges in the graph).  The row indices in each column need
00041  * not be in any particular order.  If an input column permutation is given,
00042  * node j (in the permuted matrix A*Q) is located in
00043  * Ai [Ap[Q[j]] ... Ap[Q[j]+1]-1].  This Q can be the same as the Match array
00044  * output from the maxtrans routine, for a square matrix that is structurally
00045  * full rank.
00046  *
00047  * The algorithm is from the paper by Robert E. Tarjan, "Depth-first search and
00048  * linear graph algorithms," SIAM Journal on Computing, vol. 1, no. 2,
00049  * pp. 146-160, 1972.  The time taken by strongcomp is O(nnz(A)).
00050  *
00051  * See also MC13A/B in the Harwell subroutine library (Iain S. Duff and John
00052  * K. Reid, "Algorithm 529: permutations to block triangular form," ACM Trans.
00053  * on Mathematical Software, vol. 4, no. 2, pp. 189-192, 1978, and "An
00054  * implementation of Tarjan's algorithm for the block triangular form of a
00055  * matrix," same journal, pp. 137-147.  This code is implements the same
00056  * algorithm as MC13A/B, except that the data structures are very different.
00057  * Also, unlike MC13A/B, the output permutation preserves the natural ordering
00058  * within each block.
00059  */
00060 
00061 static void amesos_dfs
00062 (
00063     /* inputs, not modified on output: */
00064     Int j,    /* start the DFS at node j */
00065     Int Ap [ ],   /* size n+1, column pointers for the matrix A */
00066     Int Ai [ ],   /* row indices, size nz = Ap [n] */
00067     Int Q [ ],    /* input column permutation */
00068 
00069     /* inputs, modified on output (each array is of size n): */
00070     Int Time [ ], /* Time [j] = "time" that node j was first visited */
00071     Int Flag [ ], /* Flag [j]: see above */
00072     Int Low [ ],  /* Low [j]: see definition below */
00073     Int *p_nblocks, /* number of blocks (aka strongly-connected-comp.)*/
00074     Int *p_timestamp, /* current "time" */
00075 
00076     /* workspace, not defined on input or output: */
00077     Int Cstack [ ], /* size n, output stack to hold nodes of components */
00078     Int Jstack [ ], /* size n, stack for the variable j */
00079     Int Pstack [ ]  /* size n, stack for the variable p */
00080 )
00081 {
00082     /* ---------------------------------------------------------------------- */
00083     /* local variables, and initializations */
00084     /* ---------------------------------------------------------------------- */
00085 
00086     /* local variables, but "global" to all DFS levels: */
00087     Int chead ;     /* top of Cstack */
00088     Int jhead ;     /* top of Jstack and Pstack */
00089 
00090     /* variables that are purely local to any one DFS level: */
00091     Int i ;     /* edge (j,i) considered; i can be next node to traverse */
00092     Int parent ;    /* parent of node j in the DFS tree */
00093     Int pend ;      /* one past the end of the adjacency list for node j */
00094     Int jj ;      /* column j of A*Q is column jj of the input matrix A */
00095 
00096     /* variables that need to be pushed then popped from the stack: */
00097     Int p ;     /* current index into the adj. list for node j */
00098     /* the variables j and p are stacked in Jstack and Pstack */
00099 
00100     /* local copies of variables in the calling routine */
00101     Int nblocks   = *p_nblocks ;
00102     Int timestamp = *p_timestamp ;
00103 
00104     /* ---------------------------------------------------------------------- */
00105     /* start a DFS at node j (same as the recursive call dfs (EMPTY, j)) */
00106     /* ---------------------------------------------------------------------- */
00107 
00108     chead = 0 ;       /* component stack is empty */
00109     jhead = 0 ;       /* Jstack and Pstack are empty */
00110     Jstack [0] = j ;      /* put the first node j on the Jstack */
00111     ASSERT (Flag [j] == UNVISITED) ;
00112 
00113     while (jhead >= 0)
00114     {
00115   j = Jstack [jhead] ;      /* grab the node j from the top of Jstack */
00116 
00117   /* determine which column jj of the A is column j of A*Q */
00118   jj = (Q == (Int *) NULL) ? (j) : (BTF_UNFLIP (Q [j])) ;
00119   pend = Ap [jj+1] ;      /* j's row index list ends at Ai [pend-1] */
00120 
00121   if (Flag [j] == UNVISITED)
00122   {
00123 
00124       /* -------------------------------------------------------------- */
00125       /* prework at node j */
00126       /* -------------------------------------------------------------- */
00127 
00128       /* node j is being visited for the first time */
00129       Cstack [++chead] = j ;      /* push j onto the stack */
00130       timestamp++ ;       /* get a timestamp */
00131       Time [j] = timestamp ;      /* give the timestamp to node j */
00132       Low [j] = timestamp ;
00133       Flag [j] = UNASSIGNED ;     /* flag node j as visited */
00134 
00135       /* -------------------------------------------------------------- */
00136       /* set Pstack [jhead] to the first entry in column j to scan */
00137       /* -------------------------------------------------------------- */
00138 
00139       Pstack [jhead] = Ap [jj] ;
00140   }
00141 
00142   /* ------------------------------------------------------------------ */
00143   /* DFS rooted at node j (start it, or continue where left off) */
00144   /* ------------------------------------------------------------------ */
00145 
00146   for (p = Pstack [jhead] ; p < pend ; p++)
00147   {
00148       i = Ai [p] ;    /* examine the edge from node j to node i */
00149       if (Flag [i] == UNVISITED)
00150       {
00151     /* Node i has not been visited - start a DFS at node i.
00152      * Keep track of where we left off in the scan of adjacency list
00153      * of node j so we can restart j where we left off. */
00154     Pstack [jhead] = p + 1 ;
00155     /* Push i onto the stack and immediately break
00156      * so we can recurse on node i. */
00157     Jstack [++jhead] = i ;
00158     ASSERT (Time [i] == EMPTY) ;
00159     ASSERT (Low [i] == EMPTY) ;
00160     /* break here to do what the recursive call dfs (j,i) does */
00161     break ;
00162       }
00163       else if (Flag [i] == UNASSIGNED)
00164       {
00165     /* Node i has been visited, but still unassigned to a block
00166      * this is a back or cross edge if Time [i] < Time [j].
00167      * Note that i might equal j, in which case this code does
00168      * nothing. */
00169     ASSERT (Time [i] > 0) ;
00170     ASSERT (Low [i] > 0) ;
00171     Low [j] = MIN (Low [j], Time [i]) ;
00172       }
00173   }
00174 
00175   if (p == pend)
00176   {
00177       /* If all adjacent nodes of j are already visited, pop j from
00178        * Jstack and do the post work for node j.  This also pops p
00179        * from the Pstack. */
00180       jhead-- ;
00181 
00182       /* -------------------------------------------------------------- */
00183       /* postwork at node j */
00184       /* -------------------------------------------------------------- */
00185 
00186       /* determine if node j is the head of a component */
00187       if (Low [j] == Time [j])
00188       {
00189     /* pop all nodes in this SCC from Cstack */
00190     while (TRUE)
00191     {
00192         ASSERT (chead >= 0) ; /* stack not empty (j in it) */
00193         i = Cstack [chead--] ;  /* pop a node from the Cstack */
00194         ASSERT (i >= 0) ;
00195         ASSERT (Flag [i] == UNASSIGNED) ;
00196         Flag [i] = nblocks ;  /* assign i to current block */
00197         if (i == j) break ;   /* current block ends at j */
00198     }
00199     nblocks++ ; /* one more block has been found */
00200       }
00201       /* update Low [parent], if the parent exists */
00202       if (jhead >= 0)
00203       {
00204     parent = Jstack [jhead] ;
00205     Low [parent] = MIN (Low [parent], Low [j]) ;
00206       }
00207   }
00208     }
00209 
00210     /* ---------------------------------------------------------------------- */
00211     /* cleanup: update timestamp and nblocks */
00212     /* ---------------------------------------------------------------------- */
00213 
00214     *p_timestamp = timestamp ;
00215     *p_nblocks   = nblocks ;
00216 }
00217 
00218 #else
00219 
00220 /* ========================================================================== */
00221 /* === dfs: recursive version (only for illustration) ======================= */
00222 /* ========================================================================== */
00223 
00224 /* The following is a recursive version of dfs, which computes identical results
00225  * as the non-recursive dfs.  It is included here because it is easier to read.
00226  * Compare the comments in the code below with the identical comments in the
00227  * non-recursive code above, and that will help you see the correlation between
00228  * the two routines.
00229  *
00230  * This routine can cause stack overflow, and is thus not recommended for heavy
00231  * usage, particularly for large matrices.  To help in delaying stack overflow,
00232  * global variables are used, reducing the amount of information each call to
00233  * dfs places on the call/return stack (the integers i, j, p, parent, and the
00234  * return address).  Note that this means the recursive code is not thread-safe.
00235  * To try this version, compile the code with -DRECURSIVE or include the
00236  * following line at the top of this file:
00237 
00238 #define RECURSIVE
00239 
00240  */
00241 
00242 static Int  /* for recursive illustration only, not for production use */
00243     chead, timestamp, nblocks, n, *Ap, *Ai, *Flag, *Cstack, *Time, *Low,
00244     *P, *R, *Q ;
00245 
00246 static void amesos_dfs
00247 (
00248     Int parent,   /* came from parent node */
00249     Int j   /* at node j in the DFS */
00250 )
00251 {
00252     Int p ;     /* current index into the adj. list for node j */
00253     Int i ;     /* edge (j,i) considered; i can be next node to traverse */
00254     Int jj ;      /* column j of A*Q is column jj of the input matrix A */
00255 
00256     /* ---------------------------------------------------------------------- */
00257     /* prework at node j */
00258     /* ---------------------------------------------------------------------- */
00259 
00260     /* node j is being visited for the first time */
00261     Cstack [++chead] = j ;      /* push j onto the stack */
00262     timestamp++ ;       /* get a timestamp */
00263     Time [j] = timestamp ;      /* give the timestamp to node j */
00264     Low [j] = timestamp ;
00265     Flag [j] = UNASSIGNED ;     /* flag node j as visited */
00266 
00267     /* ---------------------------------------------------------------------- */
00268     /* DFS rooted at node j */
00269     /* ---------------------------------------------------------------------- */
00270 
00271     /* determine which column jj of the A is column j of A*Q */
00272     jj = (Q == (Int *) NULL) ? (j) : (BTF_UNFLIP (Q [j])) ;
00273     for (p = Ap [jj] ; p < Ap [jj+1] ; p++)
00274     {
00275   i = Ai [p] ;    /* examine the edge from node j to node i */
00276   if (Flag [i] == UNVISITED)
00277   {
00278       /* Node i has not been visited - start a DFS at node i. */
00279       amesos_dfs (j, i) ;
00280   }
00281   else if (Flag [i] == UNASSIGNED)
00282   {
00283       /* Node i has been visited, but still unassigned to a block
00284        * this is a back or cross edge if Time [i] < Time [j].
00285        * Note that i might equal j, in which case this code does
00286        * nothing. */
00287       Low [j] = MIN (Low [j], Time [i]) ;
00288   }
00289     }
00290 
00291     /* ---------------------------------------------------------------------- */
00292     /* postwork at node j */
00293     /* ---------------------------------------------------------------------- */
00294 
00295     /* determine if node j is the head of a component */
00296     if (Low [j] == Time [j])
00297     {
00298   /* pop all nodes in this strongly connected component from Cstack */
00299   while (TRUE)
00300   {
00301       i = Cstack [chead--] ;  /* pop a node from the Cstack */
00302       Flag [i] = nblocks ;  /* assign node i to current block */
00303       if (i == j) break ;   /* current block ends at node j */
00304   }
00305   nblocks++ ; /* one more block has been found */
00306     }
00307     /* update Low [parent] */
00308     if (parent != EMPTY)
00309     {
00310   /* Note that this could be done with Low[j] = MIN(Low[j],Low[i]) just
00311    * after the dfs (j,i) statement above, and then parent would not have
00312    * to be an input argument.  Putting it here places all the postwork
00313    * for node j in one place, thus making the non-recursive DFS easier. */
00314   Low [parent] = MIN (Low [parent], Low [j]) ;
00315     }
00316 }
00317 
00318 #endif
00319 
00320 /* ========================================================================== */
00321 /* === btf_strongcomp ======================================================= */
00322 /* ========================================================================== */
00323 
00324 #ifndef RECURSIVE
00325 
00326 Int BTF(strongcomp) /* return # of strongly connected components */
00327 (
00328     /* input, not modified: */
00329     Int n,      /* A is n-by-n in compressed column form */
00330     Int Ap [ ],     /* size n+1 */
00331     Int Ai [ ],     /* size nz = Ap [n] */
00332 
00333     /* optional input, modified (if present) on output: */
00334     Int Q [ ],      /* size n, input column permutation.  The permutation Q can
00335          * include a flag which indicates an unmatched row.
00336          * jold = BTF_UNFLIP (Q [jnew]) is the permutation;
00337          * this function ingnores these flags.  On output, it is
00338          * modified according to the permutation P. */
00339 
00340     /* output, not defined on input: */
00341     Int P [ ],      /* size n.  P [k] = j if row and column j are kth row/col
00342          * in permuted matrix. */
00343     Int R [ ],      /* size n+1.  kth block is in rows/cols R[k] ... R[k+1]-1
00344          * of the permuted matrix. */
00345 
00346     /* workspace, not defined on input or output: */
00347     Int Work [ ]    /* size 4n */
00348 )
00349 
00350 #else
00351 
00352 Int BTF(strongcomp) /* recursive version - same as above except for Work size */
00353 (
00354     Int n_in,
00355     Int Ap_in [ ],
00356     Int Ai_in [ ],
00357     Int Q_in [ ],
00358     Int P_in [ ],
00359     Int R_in [ ],
00360     Int Work [ ]    /* size 2n */
00361 )
00362 
00363 #endif
00364 
00365 {
00366     Int j, k, b ;
00367 
00368 #ifndef RECURSIVE
00369     Int timestamp, nblocks, *Flag, *Cstack, *Time, *Low, *Jstack, *Pstack ;
00370 #else
00371     n = n_in ;
00372     Ap = Ap_in ;
00373     Ai = Ai_in ;
00374     Q = Q_in ;
00375     P = P_in ;
00376     R = R_in ;
00377     chead = EMPTY ;
00378 #endif
00379 
00380     /* ---------------------------------------------------------------------- */
00381     /* get and initialize workspace */
00382     /* ---------------------------------------------------------------------- */
00383 
00384     /* timestamp is incremented each time a new node is visited.
00385      *
00386      * Time [j] is the timestamp given to node j.
00387      *
00388      * Low [j] is the lowest timestamp of any node reachable from j via either
00389      * a path to any descendent of j in the DFS tree, or via a single edge to
00390      * an either an ancestor (a back edge) or another node that's neither an
00391      * ancestor nor a descendant (a cross edge).  If Low [j] is equal to
00392      * the timestamp of node j (Time [j]), then node j is the "head" of a
00393      * strongly connected component (SCC).  That is, it is the first node
00394      * visited in its strongly connected component, and the DFS subtree rooted
00395      * at node j spans all the nodes of the strongly connected component.
00396      *
00397      * The term "block" and "component" are used interchangebly in this code;
00398      * "block" being a matrix term and "component" being a graph term for the
00399      * same thing.
00400      *
00401      * When a node is visited, it is placed on the Cstack (for "component"
00402      * stack).  When node j is found to be an SCC head, all the nodes from the
00403      * top of the stack to node j itself form the nodes in the SCC.  This Cstack
00404      * is used for both the recursive and non-recursive versions.
00405      */
00406 
00407     Time   = Work ; Work += n ;
00408     Flag   = Work ; Work += n ;
00409     Low    = P ;    /* use output array P as workspace for Low */
00410     Cstack = R ;    /* use output array R as workspace for Cstack */
00411 
00412 #ifndef RECURSIVE
00413     /* stack for non-recursive dfs */
00414     Jstack = Work ; Work += n ;     /* stack for j */
00415     Pstack = Work ;       /* stack for p */
00416 #endif
00417 
00418     for (j = 0 ; j < n ; j++)
00419     {
00420   Flag [j] = UNVISITED ;
00421   Low [j] = EMPTY ;
00422   Time [j] = EMPTY ;
00423 #ifndef NDEBUG
00424   Cstack [j] = EMPTY ;
00425 #ifndef RECURSIVE
00426   Jstack [j] = EMPTY ;
00427   Pstack [j] = EMPTY ;
00428 #endif
00429 #endif
00430     }
00431 
00432     timestamp = 0 ; /* each node given a timestamp when it is visited */
00433     nblocks = 0 ; /* number of blocks found so far */
00434 
00435     /* ---------------------------------------------------------------------- */
00436     /* find the connected components via a depth-first-search */
00437     /* ---------------------------------------------------------------------- */
00438 
00439     for (j = 0 ; j < n ; j++)
00440     {
00441   /* node j is unvisited or assigned to a block. Cstack is empty. */
00442   ASSERT (Flag [j] == UNVISITED || (Flag [j] >= 0 && Flag [j] < nblocks));
00443   if (Flag [j] == UNVISITED)
00444   {
00445 #ifndef RECURSIVE
00446       /* non-recursive dfs (default) */
00447       amesos_dfs (j, Ap, Ai, Q, Time, Flag, Low, &nblocks, &timestamp,
00448         Cstack, Jstack, Pstack) ;
00449 #else
00450       /* recursive dfs (for illustration only) */
00451       ASSERT (chead == EMPTY) ;
00452       amesos_dfs (EMPTY, j) ;
00453       ASSERT (chead == EMPTY) ;
00454 #endif
00455   }
00456     }
00457     ASSERT (timestamp == n) ;
00458 
00459     /* ---------------------------------------------------------------------- */
00460     /* construct the block boundary array, R */
00461     /* ---------------------------------------------------------------------- */
00462 
00463     for (b = 0 ; b < nblocks ; b++)
00464     {
00465   R [b] = 0 ;
00466     }
00467     for (j = 0 ; j < n ; j++)
00468     {
00469   /* node j has been assigned to block b = Flag [j] */
00470   ASSERT (Time [j] > 0 && Time [j] <= n) ;
00471   ASSERT (Low [j] > 0 && Low [j] <= n) ;
00472   ASSERT (Flag [j] >= 0 && Flag [j] < nblocks) ;
00473   R [Flag [j]]++ ;
00474     }
00475     /* R [b] is now the number of nodes in block b.  Compute cumulative sum
00476      * of R, using Time [0 ... nblocks-1] as workspace. */
00477     Time [0] = 0 ;
00478     for (b = 1 ; b < nblocks ; b++)
00479     {
00480   Time [b] = Time [b-1] + R [b-1] ;
00481     }
00482     for (b = 0 ; b < nblocks ; b++)
00483     {
00484   R [b] = Time [b] ;
00485     }
00486     R [nblocks] = n ;
00487 
00488     /* ---------------------------------------------------------------------- */
00489     /* construct the permutation, preserving the natural order */
00490     /* ---------------------------------------------------------------------- */
00491 
00492 #ifndef NDEBUG
00493     for (k = 0 ; k < n ; k++)
00494     {
00495   P [k] = EMPTY ;
00496     }
00497 #endif
00498 
00499     for (j = 0 ; j < n ; j++)
00500     {
00501   /* place column j in the permutation */
00502   P [Time [Flag [j]]++] = j ;
00503     }
00504 
00505 #ifndef NDEBUG
00506     for (k = 0 ; k < n ; k++)
00507     {
00508   ASSERT (P [k] != EMPTY) ;
00509     }
00510 #endif
00511 
00512     /* Now block b consists of the nodes k1 to k2-1 in the permuted matrix,
00513      * where k1 = R [b] and k2 = R [b+1].  Row and column j of the original
00514      * matrix becomes row and column P [k] of the permuted matrix.  The set of
00515      * of rows/columns (nodes) in block b is given by P [k1 ... k2-1], and this
00516      * set is sorted in ascending order.  Thus, if the matrix consists of just
00517      * one block, P is the identity permutation. */
00518 
00519     /* ---------------------------------------------------------------------- */
00520     /* if Q is present on input, set Q = Q*P' */
00521     /* ---------------------------------------------------------------------- */
00522 
00523     if (Q != (Int *) NULL)
00524     {
00525   /* We found a symmetric permutation P for the matrix A*Q.  The overall
00526    * permutation is thus P*(A*Q)*P'.  Set Q=Q*P' so that the final
00527    * permutation is P*A*Q.  Use Time as workspace.  Note that this
00528    * preserves the negative values of Q if the matrix is structurally
00529    * singular. */
00530   for (k = 0 ; k < n ; k++)
00531   {
00532       Time [k] = Q [P [k]] ;
00533   }
00534   for (k = 0 ; k < n ; k++)
00535   {
00536       Q [k] = Time [k] ;
00537   }
00538     }
00539 
00540     /* ---------------------------------------------------------------------- */
00541     /* how to traverse the permuted matrix */
00542     /* ---------------------------------------------------------------------- */
00543 
00544     /* If Q is not present, the following code can be used to traverse the
00545      * permuted matrix P*A*P'
00546      *
00547      *      // compute the inverse of P
00548      *      for (knew = 0 ; knew < n ; knew++)
00549      *      {
00550      *    // row and column kold in the old matrix is row/column knew
00551      *    // in the permuted matrix P*A*P'
00552      *    kold = P [knew] ;
00553      *    Pinv [kold] = knew ;
00554      *      }
00555      *      for (b = 0 ; b < nblocks ; b++)
00556      *      {
00557      *    // traverse block b of the permuted matrix P*A*P'
00558      *    k1 = R [b] ;
00559      *    k2 = R [b+1] ;
00560      *    nk = k2 - k1 ;
00561      *    for (jnew = k1 ; jnew < k2 ; jnew++)
00562      *    {
00563      *        jold = P [jnew] ;
00564      *        for (p = Ap [jold] ; p < Ap [jold+1] ; p++)
00565      *        {
00566      *      iold = Ai [p] ;
00567      *      inew = Pinv [iold] ;
00568      *      // Entry in the old matrix is A (iold, jold), and its
00569      *      // position in the new matrix P*A*P' is (inew, jnew).
00570      *      // Let B be the bth diagonal block of the permuted
00571      *      // matrix.  If inew >= k1, then this entry is in row/
00572      *      // column (inew-k1, jnew-k1) of the nk-by-nk matrix B.
00573      *      // Otherwise, the entry is in the upper block triangular
00574      *      // part, not in any diagonal block.
00575      *        }
00576      *    }
00577      *      }
00578      *
00579      * If Q is present replace the above statement
00580      *    jold = P [jnew] ;
00581      * with
00582      *    jold = Q [jnew] ;
00583      * or
00584      *    jold = BTF_UNFLIP (Q [jnew]) ;
00585      *
00586      * then entry A (iold,jold) in the old (unpermuted) matrix is at (inew,jnew)
00587      * in the permuted matrix P*A*Q.  Everything else remains the same as the
00588      * above (simply replace P*A*P' with P*A*Q in the above comments).
00589      */
00590 
00591     /* ---------------------------------------------------------------------- */
00592     /* return # of blocks / # of strongly connected components */
00593     /* ---------------------------------------------------------------------- */
00594 
00595     return (nblocks) ;
00596 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines