Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_nesdis.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Partition/cholmod_nesdis ============================================= */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Partition Module.
00007  * Copyright (C) 2005-2006, Univ. of Florida.  Author: Timothy A. Davis
00008  * The CHOLMOD/Partition Module is licensed under Version 2.1 of the GNU
00009  * Lesser General Public License.  See lesser.txt for a text of the license.
00010  * CHOLMOD is also available under other licenses; contact authors for details.
00011  * http://www.cise.ufl.edu/research/sparse
00012  * -------------------------------------------------------------------------- */
00013 
00014 /* CHOLMOD nested dissection and graph partitioning.
00015  *
00016  * cholmod_bisect:
00017  *
00018  *  Finds a set of nodes that partitions the graph into two parts.
00019  *  Compresses the graph first.  Requires METIS.
00020  *
00021  * cholmod_nested_dissection:
00022  *
00023  *  Nested dissection, using its own compression and connected-commponents
00024  *  algorithms, an external graph partitioner (METIS), and a constrained
00025  *  minimum degree ordering algorithm (CCOLAMD or CSYMAMD).  Typically
00026  *  gives better orderings than METIS_NodeND (about 5% to 10% fewer
00027  *  nonzeros in L).
00028  *
00029  * cholmod_collapse_septree:
00030  *
00031  *  Prune the separator tree returned by cholmod_nested_dissection.
00032  *
00033  * This file contains several routines private to this file:
00034  *
00035  *  partition compress and partition a graph
00036  *  clear_flag  clear Common->Flag, but do not modify negative entries
00037  *  find_components find the connected components of a graph
00038  *
00039  * Supports any xtype (pattern, real, complex, or zomplex).
00040  */
00041 
00042 #ifndef NPARTITION
00043 
00044 #include "amesos_cholmod_internal.h"
00045 #include "amesos_cholmod_partition.h"
00046 #include "amesos_cholmod_cholesky.h"
00047 
00048 /* ========================================================================== */
00049 /* === partition ============================================================ */
00050 /* ========================================================================== */
00051 
00052 /* Find a set of nodes that partition a graph.  The graph must be symmetric
00053  * with no diagonal entries.  To compress the graph first, compress is TRUE
00054  * and on input Hash [j] holds the hash key for node j, which must be in the
00055  * range 0 to csize-1. The input graph (Cp, Ci) is destroyed.  Cew is all 1's
00056  * on input and output.  Cnw [j] > 0 is the initial weight of node j.  On
00057  * output, Cnw [i] = 0 if node i is absorbed into j and the original weight
00058  * Cnw [i] is added to Cnw [j].  If compress is FALSE, the graph is not
00059  * compressed and Cnw and Hash are unmodified.  The partition itself is held in
00060  * the output array Part of size n.  Part [j] is 0, 1, or 2, depending on
00061  * whether node j is in the left part of the graph, the right part, or the
00062  * separator, respectively.  Note that the input graph need not be connected,
00063  * and the output subgraphs (the three parts) may also be unconnected.
00064  *
00065  * Returns the size of the separator, in terms of the sum of the weights of
00066  * the nodes.  It is guaranteed to be between 1 and the total weight of all
00067  * the nodes.  If it is of size less than the total weight, then both the left
00068  * and right parts are guaranteed to be non-empty (this guarantee depends on
00069  * cholmod_metis_bisector).
00070  */
00071 
00072 static UF_long partition    /* size of separator or -1 if failure */
00073 (
00074     /* inputs, not modified on output */
00075 #ifndef NDEBUG
00076     Int csize,    /* upper bound on # of edges in the graph;
00077        * csize >= MAX (n, nnz(C)) must hold. */
00078 #endif
00079     int compress, /* if TRUE the compress the graph first */
00080 
00081     /* input/output */
00082     Int Hash [ ], /* Hash [i] = hash >= 0 is the hash function for node
00083        * i on input.  On output, Hash [i] = FLIP (j) if node
00084        * i is absorbed into j.  Hash [i] >= 0 if i has not
00085        * been absorbed. */
00086 
00087     /* input graph, compressed graph of cn nodes on output */
00088     cholmod_sparse *C,
00089 
00090     /* input/output */
00091     Int Cnw [ ],  /* size n.  Cnw [j] > 0 is the weight of node j on
00092        * input.  On output, if node i is absorbed into
00093        * node j, then Cnw [i] = 0 and the original weight of
00094        * node i is added to Cnw [j].  The sum of Cnw [0..n-1]
00095        * is not modified. */
00096 
00097     /* workspace */
00098     Int Cew [ ],  /* size csize, all 1's on input and output */
00099 
00100     /* more workspace, undefined on input and output */
00101     Int Cmap [ ], /* size n (i/i/l) */
00102 
00103     /* output */
00104     Int Part [ ], /* size n, Part [j] = 0, 1, or 2. */
00105 
00106     cholmod_common *Common
00107 )
00108 {
00109     Int n, hash, head, i, j, k, p, pend, ilen, ilast, pi, piend,
00110   jlen, ok, cn, csep, pdest, nodes_pruned, nz, total_weight, jscattered ;
00111     Int *Cp, *Ci, *Next, *Hhead ;
00112 
00113 #ifndef NDEBUG
00114     Int cnt, pruned ;
00115     double work = 0, goodwork = 0 ;
00116 #endif
00117 
00118     /* ---------------------------------------------------------------------- */
00119     /* quick return for small or empty graphs */
00120     /* ---------------------------------------------------------------------- */
00121 
00122     n = C->nrow ;
00123     Cp = C->p ;
00124     Ci = C->i ;
00125     nz = Cp [n] ;
00126 
00127     PRINT2 (("Partition start, n "ID" nz "ID"\n", n, nz)) ;
00128 
00129     total_weight = 0 ;
00130     for (j = 0 ; j < n ; j++)
00131     {
00132   ASSERT (Cnw [j] > 0) ;
00133   total_weight += Cnw [j] ;
00134     }
00135 
00136     if (n <= 2)
00137     {
00138   /* very small graph */
00139   for (j = 0 ; j < n ; j++)
00140   {
00141       Part [j] = 2 ;
00142   }
00143   return (total_weight) ;
00144     }
00145     else if (nz <= 0)
00146     {
00147   /* no edges, this is easy */
00148   PRINT2 (("diagonal matrix\n")) ;
00149   k = n/2 ;
00150   for (j = 0 ; j < k ; j++)
00151   {
00152       Part [j] = 0 ;
00153   }
00154   for ( ; j < n ; j++)
00155   {
00156       Part [j] = 1 ;
00157   }
00158   /* ensure the separator is not empty (required by nested dissection) */
00159   Part [n-1] = 2 ;
00160   return (Cnw [n-1]) ;
00161     }
00162 
00163 #ifndef NDEBUG
00164     ASSERT (n > 1 && nz > 0) ;
00165     PRINT2 (("original graph:\n")) ;
00166     for (j = 0 ; j < n ; j++)
00167     {
00168   PRINT2 ((""ID": ", j)) ;
00169   for (p = Cp [j] ; p < Cp [j+1] ; p++)
00170   {
00171       i = Ci [p] ;
00172       PRINT3 ((""ID" ", i)) ;
00173       ASSERT (i >= 0 && i < n && i != j) ;
00174   }
00175   PRINT2 (("hash: "ID"\n", Hash [j])) ;
00176     }
00177     DEBUG (for (p = 0 ; p < csize ; p++) ASSERT (Cew [p] == 1)) ;
00178 #endif
00179 
00180     nodes_pruned = 0 ;
00181 
00182     if (compress)
00183     {
00184 
00185   /* ------------------------------------------------------------------ */
00186   /* get workspace */
00187   /* ------------------------------------------------------------------ */
00188 
00189   Next = Part ; /* use Part as workspace for Next [ */
00190   Hhead = Cew ; /* use Cew as workspace for Hhead [ */
00191 
00192   /* ------------------------------------------------------------------ */
00193   /* create the hash buckets */
00194   /* ------------------------------------------------------------------ */
00195 
00196   for (j = 0 ; j < n ; j++)
00197   {
00198       /* get the hash key for node j */
00199       hash = Hash [j] ;
00200       ASSERT (hash >= 0 && hash < csize) ;
00201       head = Hhead [hash] ;
00202       if (head > EMPTY)
00203       {
00204     /* hash bucket for this hash key is empty. */
00205     head = EMPTY ;
00206       }
00207       else
00208       {
00209     /* hash bucket for this hash key is not empty.  get old head */
00210     head = FLIP (head) ;
00211     ASSERT (head >= 0 && head < n) ;
00212       }
00213       /* node j becomes the new head of the hash bucket.  FLIP it so that
00214        * we can tell the difference between an empty or non-empty hash
00215        * bucket. */
00216       Hhead [hash] = FLIP (j) ;
00217       Next [j] = head ;
00218       ASSERT (head >= EMPTY && head < n) ;
00219   }
00220 
00221 #ifndef NDEBUG
00222   for (cnt = 0, k = 0 ; k < n ; k++)
00223   {
00224       ASSERT (Hash [k] >= 0 && Hash [k] < csize) ;    /* k is alive */
00225       hash = Hash [k] ;
00226       ASSERT (hash >= 0 && hash < csize) ;
00227       head = Hhead [hash] ;
00228       ASSERT (head < EMPTY) ; /* hash bucket not empty */
00229       j = FLIP (head) ;
00230       ASSERT (j >= 0 && j < n) ;
00231       if (j == k)
00232       {
00233     PRINT2 (("hash "ID": ", hash)) ;
00234     for ( ; j != EMPTY ; j = Next [j])
00235     {
00236         PRINT3 ((" "ID"", j)) ;
00237         ASSERT (j >= 0 && j < n) ;
00238         ASSERT (Hash [j] == hash) ;
00239         cnt++ ;
00240         ASSERT (cnt <= n) ;
00241     }
00242     PRINT2 (("\n")) ;
00243       }
00244   }
00245   ASSERT (cnt == n) ;
00246 #endif
00247 
00248   /* ------------------------------------------------------------------ */
00249   /* scan the non-empty hash buckets for indistinguishable nodes */
00250   /* ------------------------------------------------------------------ */
00251 
00252   /* If there are no hash collisions and no compression occurs, this takes
00253    * O(n) time.  If no hash collisions, but some nodes are removed, this
00254    * takes time O(n+e) where e is the sum of the degress of the nodes
00255    * that are removed.  Even with many hash collisions (a rare case),
00256    * this algorithm has never been observed to perform more than nnz(A)
00257    * useless work.
00258    *
00259    * Cmap is used as workspace to mark nodes of the graph, [
00260    * for comparing the nonzero patterns of two nodes i and j.
00261    */
00262 
00263 #define Cmap_MARK(i)   Cmap [i] = j
00264 #define Cmap_MARKED(i) (Cmap [i] == j)
00265 
00266   for (i = 0 ; i < n ; i++)
00267   {
00268       Cmap [i] = EMPTY ;
00269   }
00270 
00271   for (k = 0 ; k < n ; k++)
00272   {
00273       hash = Hash [k] ;
00274       ASSERT (hash >= FLIP (n-1) && hash < csize) ;
00275       if (hash < 0)
00276       {
00277     /* node k has already been absorbed into some other node */
00278     ASSERT (FLIP (Hash [k]) >= 0 && FLIP (Hash [k] < n)) ;
00279     continue ;
00280       }
00281       head = Hhead [hash] ;
00282       ASSERT (head < EMPTY || head == 1) ;
00283       if (head == 1)
00284       {
00285     /* hash bucket is already empty */
00286     continue ;
00287       }
00288       PRINT2 (("\n--------------------hash "ID":\n", hash)) ;
00289       for (j = FLIP (head) ; j != EMPTY && Next[j] > EMPTY ; j = Next [j])
00290       {
00291     /* compare j with all nodes i following it in hash bucket */
00292     ASSERT (j >= 0 && j < n && Hash [j] == hash) ;
00293     p = Cp [j] ;
00294     pend = Cp [j+1] ;
00295     jlen = pend - p ;
00296     jscattered = FALSE ;
00297     DEBUG (for (i = 0 ; i < n ; i++) ASSERT (!Cmap_MARKED (i))) ;
00298     DEBUG (pruned = FALSE) ;
00299     ilast = j ;
00300     for (i = Next [j] ; i != EMPTY ; i = Next [i])
00301     {
00302         ASSERT (i >= 0 && i < n && Hash [i] == hash && i != j) ;
00303         pi = Cp [i] ;
00304         piend = Cp [i+1] ;
00305         ilen = piend - pi ;
00306         DEBUG (work++) ;
00307         if (ilen != jlen)
00308         {
00309       /* i and j have different degrees */
00310       ilast = i ;
00311       continue ;
00312         }
00313         /* scatter the pattern of node j, if not already */
00314         if (!jscattered)
00315         {
00316       Cmap_MARK (j) ;
00317       for ( ; p < pend ; p++)
00318       {
00319           Cmap_MARK (Ci [p]) ;
00320       }
00321       jscattered = TRUE ;
00322       DEBUG (work += jlen) ;
00323         }
00324         for (ok = Cmap_MARKED (i) ; ok && pi < piend ; pi++)
00325         {
00326       ok = Cmap_MARKED (Ci [pi]) ;
00327       DEBUG (work++) ;
00328         }
00329         if (ok)
00330         {
00331       /* found it.  kill node i and merge it into j */
00332       PRINT2 (("found "ID" absorbed into "ID"\n", i, j)) ;
00333       Hash [i] = FLIP (j) ;
00334       Cnw [j] += Cnw [i] ;
00335       Cnw [i] = 0 ;
00336       ASSERT (ilast != i && ilast >= 0 && ilast < n) ;
00337       Next [ilast] = Next [i] ; /* delete i from bucket */
00338       nodes_pruned++ ;
00339       DEBUG (goodwork += (ilen+1)) ;
00340       DEBUG (pruned = TRUE) ;
00341         }
00342         else
00343         {
00344       /* i and j are different */
00345       ilast = i ;
00346         }
00347     }
00348     DEBUG (if (pruned) goodwork += jlen) ;
00349       }
00350       /* empty the hash bucket, restoring Cew */
00351       Hhead [hash] = 1 ;
00352   }
00353 
00354   DEBUG (if (((work - goodwork) / (double) nz) > 0.20) PRINT0 ((
00355       "work %12g good %12g nz %12g (wasted work/nz: %6.2f )\n",
00356       work, goodwork, (double) nz, (work - goodwork) / ((double) nz)))) ;
00357 
00358   /* All hash buckets now empty.  Cmap no longer needed as workspace. ]
00359    * Cew no longer needed as Hhead; Cew is now restored to all ones. ]
00360    * Part no longer needed as workspace for Next. ] */
00361     }
00362 
00363     /* Edge weights are all one, node weights reflect node absorption */
00364     DEBUG (for (p = 0 ; p < csize ; p++) ASSERT (Cew [p] == 1)) ;
00365     DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) cnt += Cnw [j]) ;
00366     ASSERT (cnt == total_weight) ;
00367 
00368     /* ---------------------------------------------------------------------- */
00369     /* compress and partition the graph */
00370     /* ---------------------------------------------------------------------- */
00371 
00372     if (nodes_pruned == 0)
00373     {
00374 
00375   /* ------------------------------------------------------------------ */
00376   /* no pruning done at all.  Do not create the compressed graph */
00377   /* ------------------------------------------------------------------ */
00378 
00379   /* FUTURE WORK: could call CHACO, SCOTCH, ... here too */
00380   csep = CHOLMOD(metis_bisector) (C, Cnw, Cew, Part, Common) ;
00381 
00382     }
00383     else if (nodes_pruned == n-1)
00384     {
00385 
00386   /* ------------------------------------------------------------------ */
00387   /* only one node left.  This is a dense graph */
00388   /* ------------------------------------------------------------------ */
00389 
00390   PRINT2 (("completely dense graph\n")) ;
00391   csep = total_weight ;
00392   for (j = 0 ; j < n ; j++)
00393   {
00394       Part [j] = 2 ;
00395   }
00396 
00397     }
00398     else
00399     {
00400 
00401   /* ------------------------------------------------------------------ */
00402   /* compress the graph and partition the compressed graph */
00403   /* ------------------------------------------------------------------ */
00404 
00405   /* ------------------------------------------------------------------ */
00406   /* create the map from the uncompressed graph to the compressed graph */
00407   /* ------------------------------------------------------------------ */
00408 
00409   /* Cmap [j] = k if node j is alive and the kth node of compressed graph.
00410    * The mapping is done monotonically (that is, k <= j) to simplify the
00411    * uncompression later on.  Cmap [j] = EMPTY if node j is dead. */
00412 
00413   for (j = 0 ; j < n ; j++)
00414   {
00415       Cmap [j] = EMPTY ;
00416   }
00417   k = 0 ;
00418   for (j = 0 ; j < n ; j++)
00419   {
00420       if (Cnw [j] > 0)
00421       {
00422     ASSERT (k <= j) ;
00423     Cmap [j] = k++ ;
00424       }
00425   }
00426   cn = k ;      /* # of nodes in compressed graph */
00427   PRINT2 (("compressed graph from "ID" to "ID" nodes\n", n, cn)) ;
00428   ASSERT (cn > 1 && cn == n - nodes_pruned) ;
00429 
00430   /* ------------------------------------------------------------------ */
00431   /* create the compressed graph */
00432   /* ------------------------------------------------------------------ */
00433 
00434   k = 0 ;
00435   pdest = 0 ;
00436   for (j = 0 ; j < n ; j++)
00437   {
00438       if (Cnw [j] > 0)
00439       {
00440     /* node j in the full graph is node k in the compressed graph */
00441     ASSERT (k <= j && Cmap [j] == k) ;
00442     p = Cp [j] ;
00443     pend = Cp [j+1] ;
00444     Cp [k] = pdest ;
00445     Cnw [k] = Cnw [j] ;
00446     for ( ; p < pend ; p++)
00447     {
00448         /* prune dead nodes, and remap to new node numbering */
00449         i = Ci [p] ;
00450         ASSERT (i >= 0 && i < n && i != j) ;
00451         i = Cmap [i] ;
00452         ASSERT (i >= EMPTY && i < cn && i != k) ;
00453         if (i > EMPTY)
00454         {
00455       ASSERT (pdest <= p) ;
00456       Ci [pdest++] = i ;
00457         }
00458     }
00459     k++ ;
00460       }
00461   }
00462   Cp [cn] = pdest ;
00463   C->nrow = cn ;
00464   C->ncol = cn ;  /* affects mem stats unless restored when C free'd */
00465 
00466 #ifndef NDEBUG
00467   PRINT2 (("pruned graph ("ID"/"ID") nodes, ("ID"/"ID") edges\n",
00468         cn, n, pdest, nz)) ;
00469   PRINT2 (("compressed graph:\n")) ;
00470   for (cnt = 0, j = 0 ; j < cn ; j++)
00471   {
00472       PRINT2 ((""ID": ", j)) ;
00473       for (p = Cp [j] ; p < Cp [j+1] ; p++)
00474       {
00475     i = Ci [p] ;
00476     PRINT3 ((""ID" ", i)) ;
00477     ASSERT (i >= 0 && i < cn && i != j) ;
00478       }
00479       PRINT2 (("weight: "ID"\n", Cnw [j])) ;
00480       ASSERT (Cnw [j] > 0) ;
00481       cnt += Cnw [j] ;
00482   }
00483   ASSERT (cnt == total_weight) ;
00484   for (j = 0 ; j < n ; j++) PRINT2 (("Cmap ["ID"] = "ID"\n", j, Cmap[j]));
00485   ASSERT (k == cn) ;
00486 #endif
00487 
00488   /* ------------------------------------------------------------------ */
00489   /* find the separator of the compressed graph */
00490   /* ------------------------------------------------------------------ */
00491 
00492   /* FUTURE WORK: could call CHACO, SCOTCH, ... here too */
00493   csep = CHOLMOD(metis_bisector) (C, Cnw, Cew, Part, Common) ;
00494 
00495   if (csep < 0)
00496   {
00497       /* failed */
00498       return (-1) ;
00499   }
00500 
00501   PRINT2 (("Part: ")) ;
00502   DEBUG (for (j = 0 ; j < cn ; j++) PRINT2 ((""ID" ", Part [j]))) ;
00503   PRINT2 (("\n")) ;
00504 
00505   /* Cp and Ci no longer needed */
00506 
00507   /* ------------------------------------------------------------------ */
00508   /* find the separator of the uncompressed graph */
00509   /* ------------------------------------------------------------------ */
00510 
00511   /* expand the separator to live nodes in the uncompressed graph */
00512   for (j = n-1 ; j >= 0 ; j--)
00513   {
00514       /* do this in reverse order so that Cnw can be expanded in place */
00515       k = Cmap [j] ;
00516       ASSERT (k >= EMPTY && k < n) ;
00517       if (k > EMPTY)
00518       {
00519     /* node k in compressed graph and is node j in full graph */
00520     ASSERT (k <= j) ;
00521     ASSERT (Hash [j] >= EMPTY) ;
00522     Part [j] = Part [k] ;
00523     Cnw [j] = Cnw [k] ;
00524       }
00525       else
00526       {
00527     /* node j is a dead node */
00528     Cnw [j] = 0 ;
00529     DEBUG (Part [j] = EMPTY) ;
00530     ASSERT (Hash [j] < EMPTY) ;
00531       }
00532   }
00533 
00534   /* find the components for the dead nodes */
00535   for (i = 0 ; i < n ; i++)
00536   {
00537       if (Hash [i] < EMPTY)
00538       {
00539     /* node i has been absorbed into node j */
00540     j = FLIP (Hash [i]) ;
00541     ASSERT (Part [i] == EMPTY && j >= 0 && j < n && Cnw [i] == 0) ;
00542     Part [i] = Part [j] ;
00543       }
00544       ASSERT (Part [i] >= 0 && Part [i] <= 2) ;
00545   }
00546 
00547 #ifndef NDEBUG
00548   PRINT2 (("Part: ")) ;
00549   for (cnt = 0, j = 0 ; j < n ; j++)
00550   {
00551       ASSERT (Part [j] != EMPTY) ;
00552       PRINT2 ((""ID" ", Part [j])) ;
00553       if (Part [j] == 2) cnt += Cnw [j] ;
00554   }
00555   PRINT2 (("\n")) ;
00556   PRINT2 (("csep "ID" "ID"\n", cnt, csep)) ;
00557   ASSERT (cnt == csep) ;
00558   for (cnt = 0, j = 0 ; j < n ; j++) cnt += Cnw [j] ;
00559   ASSERT (cnt == total_weight) ;
00560 #endif
00561 
00562     }
00563 
00564     /* ---------------------------------------------------------------------- */
00565     /* return the separator (or -1 if error) */
00566     /* ---------------------------------------------------------------------- */
00567 
00568     PRINT2 (("Partition done, n "ID" csep "ID"\n", n, csep)) ;
00569     return (csep) ;
00570 }
00571 
00572 
00573 /* ========================================================================== */
00574 /* === clear_flag =========================================================== */
00575 /* ========================================================================== */
00576 
00577 /* A node j has been removed from the graph if Flag [j] < EMPTY.
00578  * If Flag [j] >= EMPTY && Flag [j] < mark, then node j is alive but unmarked.
00579  * Flag [j] == mark means that node j is alive and marked.  Incrementing mark
00580  * means that all nodes are either (still) dead, or live but unmarked.
00581  *
00582  * On output, Common->mark < Common->Flag [i] for all i from 0 to Common->nrow.
00583  * This is the same output condition as cholmod_clear_flag, except that this
00584  * routine maintains the Flag [i] < EMPTY condition as well, if that condition
00585  * was true on input.
00586  *
00587  * workspace: Flag (nrow)
00588  */
00589 
00590 static UF_long clear_flag (cholmod_common *Common)
00591 {
00592     Int nrow, i ;
00593     Int *Flag ;
00594     PRINT2 (("old mark %ld\n", Common->mark)) ;
00595     Common->mark++ ;
00596     PRINT2 (("new mark %ld\n", Common->mark)) ;
00597     if (Common->mark <= 0)
00598     {
00599   nrow = Common->nrow ;
00600   Flag = Common->Flag ;
00601   for (i = 0 ; i < nrow ; i++)
00602   {
00603       /* if Flag [i] < EMPTY, leave it alone */
00604       if (Flag [i] >= EMPTY)
00605       {
00606     Flag [i] = EMPTY ;
00607       }
00608   }
00609   /* now Flag [i] <= EMPTY for all i */
00610   Common->mark = 0 ;
00611     }
00612     return (Common->mark) ;
00613 }
00614 
00615 
00616 /* ========================================================================== */
00617 /* === find_components ====================================================== */
00618 /* ========================================================================== */
00619 
00620 /* Find all connected components of the current subgraph C.  The subgraph C
00621  * consists of the nodes of B that appear in the set Map [0..cn-1].  If Map
00622  * is NULL, then it is assumed to be the identity mapping
00623  * (Map [0..cn-1] = 0..cn-1).
00624  *
00625  * A node j does not appear in B if it has been ordered (Flag [j] < EMPTY,
00626  * which means that j has been ordered and is "deleted" from B).
00627  *
00628  * If the size of a component is large, it is placed on the component stack,
00629  * Cstack.  Otherwise, its nodes are ordered and it is not placed on the Cstack.
00630  *
00631  * A component S is defined by a "representative node" (repnode for short)
00632  * called the snode, which is one of the nodes in the subgraph.  Likewise, the
00633  * subgraph C is defined by its repnode, called cnode.
00634  * 
00635  * If Part is not NULL on input, then Part [i] determines how the components
00636  * are placed on the stack.  Components containing nodes i with Part [i] == 0
00637  * are placed first, followed by components with Part [i] == 1. 
00638  *
00639  * The first node placed in each of the two parts is flipped when placed in
00640  * the Cstack.  This allows the components of the two parts to be found simply
00641  * by traversing the Cstack.
00642  *
00643  * workspace: Flag (nrow)
00644  */
00645 
00646 static void find_components
00647 (
00648     /* inputs, not modified on output */
00649     cholmod_sparse *B,
00650     Int Map [ ],      /* size n, only Map [0..cn-1] used */
00651     Int cn,       /* # of nodes in C */
00652     Int cnode,        /* root node of component C, or EMPTY if C is the
00653            * entire graph B */
00654 
00655     Int Part [ ],     /* size cn, optional */
00656 
00657     /* input/output */
00658     Int Bnz [ ],      /* size n.  Bnz [j] = # nonzeros in column j of B.
00659            * Reduce since B is pruned of dead nodes. */
00660 
00661     Int CParent [ ],      /* CParent [i] = j if component with repnode j is
00662            * the parent of the component with repnode i.
00663            * CParent [i] = EMPTY if the component with
00664            * repnode i is a root of the separator tree.
00665            * CParent [i] is -2 if i is not a repnode. */
00666     Int Cstack [ ],     /* component stack for nested dissection */
00667     Int *top,       /* Cstack [0..top] contains root nodes of the
00668            * the components currently in the stack */
00669 
00670     /* workspace, undefined on input and output: */
00671     Int Queue [ ],      /* size n, for breadth-first search */
00672 
00673     cholmod_common *Common
00674 )
00675 {
00676     Int n, mark, cj, j, sj, sn, p, i, snode, pstart, pdest, pend, nd_components,
00677   part, first ;
00678     Int *Bp, *Bi, *Flag ;
00679 
00680     /* ---------------------------------------------------------------------- */
00681     /* get workspace */
00682     /* ---------------------------------------------------------------------- */
00683 
00684     PRINT2 (("find components: cn %d\n", cn)) ;
00685     Flag = Common->Flag ;     /* size n */
00686     Common->mark = EMPTY ;      /* force initialization of Flag array */
00687     mark = clear_flag (Common) ;    /* clear Flag but preserve Flag [i]<EMPTY */
00688     Bp = B->p ;
00689     Bi = B->i ;
00690     n = B->nrow ;
00691     ASSERT (cnode >= EMPTY && cnode < n) ;
00692     ASSERT (IMPLIES (cnode >= 0, Flag [cnode] < EMPTY)) ;
00693 
00694     /* get ordering parameters */
00695     nd_components = Common->method [Common->current].nd_components ;
00696 
00697     /* ---------------------------------------------------------------------- */
00698     /* find the connected components of C via a breadth-first search */
00699     /* ---------------------------------------------------------------------- */
00700 
00701     part = (Part == NULL) ? 0 : 1 ;
00702 
00703     /* examine each part (part 1 and then part 0) */
00704     for (part = (Part == NULL) ? 0 : 1 ; part >= 0 ; part--)
00705     {
00706 
00707   /* first is TRUE for the first connected component in each part */
00708   first = TRUE ;
00709 
00710   /* find all connected components in the current part */
00711   for (cj = 0 ; cj < cn ; cj++)
00712   {
00713       /* get node snode, which is node cj of C.  It might already be in
00714        * the separator of C (and thus ordered, with Flag [snode] < EMPTY)
00715        */
00716       snode = (Map == NULL) ? (cj) : (Map [cj]) ;
00717       ASSERT (snode >= 0 && snode < n) ;
00718 
00719       if (Flag [snode] >= EMPTY && Flag [snode] < mark
00720         && ((Part == NULL) || Part [cj] == part))
00721       {
00722 
00723     /* ---------------------------------------------------------- */
00724     /* find new connected component S */
00725     /* ---------------------------------------------------------- */
00726 
00727     /* node snode is the repnode of a connected component S, the
00728      * parent of which is cnode, the repnode of C.  If cnode is
00729      * EMPTY then C is the original graph B. */
00730     PRINT2 (("----------:::snode "ID" cnode "ID"\n", snode, cnode));
00731 
00732     ASSERT (CParent [snode] == -2) ;
00733     if (first || nd_components)
00734     {
00735         /* If this is the first node in this part, then it becomes
00736          * the repnode of all components in this part, and all
00737          * components in this part form a single node in the
00738          * separator tree.  If nd_components is TRUE, then all
00739          * connected components form their own node in the
00740          * separator tree.
00741          */
00742         CParent [snode] = cnode ;
00743     }
00744 
00745     /* place j in the queue and mark it */
00746     sj = 0 ;
00747     Queue [0] = snode ;
00748     Flag [snode] = mark ;
00749     sn = 1 ;
00750 
00751     /* breadth-first traversal, starting at node j */
00752     for (sj = 0 ; sj < sn ; sj++)
00753     {
00754         /* get node j from head of Queue and traverse its edges */
00755         j = Queue [sj] ;
00756         PRINT2 (("    j: "ID"\n", j)) ;
00757         ASSERT (j >= 0 && j < n) ;
00758         ASSERT (Flag [j] == mark) ;
00759         pstart = Bp [j] ;
00760         pdest = pstart ;
00761         pend = pstart + Bnz [j] ;
00762         for (p = pstart ; p < pend ; p++)
00763         {
00764       i = Bi [p] ;
00765       if (i != j && Flag [i] >= EMPTY)
00766       {
00767           /* node is still in the graph */
00768           Bi [pdest++] = i ;
00769           if (Flag [i] < mark)
00770           {
00771         /* node i is in this component S, and unflagged
00772          * (first time node i has been seen in this BFS).
00773          * place node i in the queue and mark it */
00774         Queue [sn++] = i ;
00775         Flag [i] = mark ;
00776           }
00777       }
00778         }
00779         /* edges to dead nodes have been removed */
00780         Bnz [j] = pdest - pstart ;
00781     }
00782 
00783     /* ---------------------------------------------------------- */
00784     /* order S if it is small; place it on Cstack otherwise */
00785     /* ---------------------------------------------------------- */
00786 
00787     PRINT2 (("sn "ID"\n", sn)) ;
00788 
00789     /* place the new component on the Cstack.  Flip the node if
00790      * is the first connected component of the current part,
00791      * or if all components are treated as their own node in
00792      * the separator tree. */
00793     Cstack [++(*top)] =
00794       (first || nd_components) ? FLIP (snode) : snode ;
00795     first = FALSE ;
00796       }
00797   }
00798     }
00799 
00800     /* clear Flag array, but preserve Flag [i] < EMPTY */
00801     clear_flag (Common) ;
00802 }
00803 
00804 
00805 /* ========================================================================== */
00806 /* === cholmod_bisect ======================================================= */
00807 /* ========================================================================== */
00808 
00809 /* Finds a node bisector of A, A*A', A(:,f)*A(:,f)'.
00810  *
00811  * workspace: Flag (nrow),
00812  *  Iwork (nrow if symmetric, max (nrow,ncol) if unsymmetric).
00813  *  Allocates a temporary matrix B=A*A' or B=A,
00814  *  and O(nnz(A)) temporary memory space.
00815  */
00816 
00817 UF_long CHOLMOD(bisect) /* returns # of nodes in separator */
00818 (
00819     /* ---- input ---- */
00820     cholmod_sparse *A,  /* matrix to bisect */
00821     Int *fset,    /* subset of 0:(A->ncol)-1 */
00822     size_t fsize, /* size of fset */
00823     int compress, /* if TRUE, compress the graph first */
00824     /* ---- output --- */
00825     Int *Partition, /* size A->nrow.  Node i is in the left graph if
00826        * Partition [i] = 0, the right graph if 1, and in the
00827        * separator if 2. */
00828     /* --------------- */
00829     cholmod_common *Common
00830 )
00831 {
00832     Int *Bp, *Bi, *Hash, *Cmap, *Bnw, *Bew, *Iwork ;
00833     cholmod_sparse *B ;
00834     unsigned Int hash ;
00835     Int j, n, bnz, sepsize, p, pend ;
00836     size_t csize, s ;
00837     int ok = TRUE ;
00838 
00839     /* ---------------------------------------------------------------------- */
00840     /* check inputs */
00841     /* ---------------------------------------------------------------------- */
00842 
00843     RETURN_IF_NULL_COMMON (EMPTY) ;
00844     RETURN_IF_NULL (A, EMPTY) ;
00845     RETURN_IF_NULL (Partition, EMPTY) ;
00846     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
00847     Common->status = CHOLMOD_OK ;
00848 
00849     /* ---------------------------------------------------------------------- */
00850     /* quick return */
00851     /* ---------------------------------------------------------------------- */
00852 
00853     n = A->nrow ;
00854     if (n == 0)
00855     {
00856   return (0) ;
00857     }
00858 
00859     /* ---------------------------------------------------------------------- */
00860     /* allocate workspace */
00861     /* ---------------------------------------------------------------------- */
00862 
00863     /* s = n + MAX (n, A->ncol) */
00864     s = CHOLMOD(add_size_t) (A->nrow, MAX (A->nrow, A->ncol), &ok) ;
00865     if (!ok)
00866     {
00867   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00868   return (EMPTY) ;
00869     }
00870 
00871     CHOLMOD(allocate_work) (n, s, 0, Common) ;
00872     if (Common->status < CHOLMOD_OK)
00873     {
00874   return (EMPTY) ;
00875     }
00876     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00877 
00878     Iwork = Common->Iwork ;
00879     Hash = Iwork ;    /* size n, (i/l/l) */
00880     Cmap = Iwork + n ;    /* size n, (i/i/l) */
00881 
00882     /* ---------------------------------------------------------------------- */
00883     /* convert the matrix to adjacency list form */
00884     /* ---------------------------------------------------------------------- */
00885 
00886     /* The input graph to must be symmetric, with no diagonal entries
00887      * present.  The columns need not be sorted. */
00888 
00889     /* B = A, A*A', or A(:,f)*A(:,f)', upper and lower parts present */
00890 
00891     if (A->stype)
00892     {
00893   /* Add the upper/lower part to a symmetric lower/upper matrix by
00894    * converting to unsymmetric mode */
00895   /* workspace: Iwork (nrow) */
00896   B = CHOLMOD(copy) (A, 0, -1, Common) ;
00897     }
00898     else
00899     {
00900   /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
00901   /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
00902   B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
00903     }
00904 
00905     if (Common->status < CHOLMOD_OK)
00906     {
00907   return (EMPTY) ;
00908     }
00909     Bp = B->p ;
00910     Bi = B->i ;
00911     bnz = Bp [n] ;
00912     ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;
00913 
00914     /* B does not include the diagonal, and both upper and lower parts.
00915      * Common->anz includes the diagonal, and just the lower part of B */
00916     Common->anz = bnz / 2 + ((double) n) ;
00917 
00918     /* Bew should be at least size n for the hash function to work well */
00919     /* this cannot cause overflow, because the matrix is already created */
00920     csize = MAX (((size_t) n) + 1, (size_t) bnz) ;
00921 
00922     /* create the graph using Flag as workspace for node weights [ */
00923     Bnw = Common->Flag ;    /* size n workspace */
00924 
00925     /* compute hash for each node if compression requested */
00926     if (compress)
00927     {
00928   for (j = 0 ; j < n ; j++)
00929   {
00930       hash = j ;
00931       pend = Bp [j+1] ;
00932       for (p = Bp [j] ; p < pend ; p++)
00933       {
00934     hash += Bi [p] ;
00935     ASSERT (Bi [p] != j) ;
00936       }
00937       /* finalize the hash key for node j */
00938       hash %= csize ;
00939       Hash [j] = (Int) hash ;
00940       ASSERT (Hash [j] >= 0 && Hash [j] < csize) ;
00941   }
00942     }
00943 
00944     /* allocate edge weights */
00945     Bew = CHOLMOD(malloc) (csize, sizeof (Int), Common) ;
00946     if (Common->status < CHOLMOD_OK)
00947     {
00948   /* out of memory */
00949   CHOLMOD(free_sparse) (&B, Common) ;
00950   CHOLMOD(free) (csize, sizeof (Int), Bew, Common) ;
00951   return (EMPTY) ;
00952     }
00953 
00954     /* graph has unit node and edge weights */
00955     for (j = 0 ; j < n ; j++)
00956     {
00957   Bnw [j] = 1 ;
00958     }
00959     for (s = 0 ; s < csize ; s++)
00960     {
00961   Bew [s] = 1 ;
00962     }
00963 
00964     /* ---------------------------------------------------------------------- */
00965     /* compress and partition the graph */
00966     /* ---------------------------------------------------------------------- */
00967 
00968     sepsize = partition (
00969 #ifndef NDEBUG
00970       csize,
00971 #endif
00972       compress, Hash, B, Bnw, Bew, Cmap, Partition, Common) ;
00973 
00974     /* contents of Bp, Bi, Bnw, and Bew no longer needed ] */
00975 
00976     /* If partition fails, free the workspace below and return sepsize < 0 */
00977 
00978     /* ---------------------------------------------------------------------- */
00979     /* free workspace */
00980     /* ---------------------------------------------------------------------- */
00981 
00982     B->ncol = n ;   /* restore size for memory usage statistics */
00983     CHOLMOD(free_sparse) (&B, Common) ;
00984     Common->mark = EMPTY ;
00985     CHOLMOD(clear_flag) (Common) ;
00986     CHOLMOD(free) (csize, sizeof (Int), Bew, Common) ;
00987     return (sepsize) ;
00988 }
00989 
00990 
00991 /* ========================================================================== */
00992 /* === cholmod_nested_dissection ============================================ */
00993 /* ========================================================================== */
00994 
00995 /* This method uses a node bisector, applied recursively (but using a
00996  * non-recursive algorithm).  Once the graph is partitioned, it calls a
00997  * constrained min degree code (CAMD or CSYMAMD for A+A', and CCOLAMD for A*A')
00998  * to order all the nodes in the graph - but obeying the constraints determined
00999  * by the separators.  This routine is similar to METIS_NodeND, except for how
01000  * it treats the leaf nodes.  METIS_NodeND orders the leaves of the separator
01001  * tree with MMD, ignoring the rest of the matrix when ordering a single leaf.
01002  * This routine orders the whole matrix with CSYMAMD or CCOLAMD, all at once,
01003  * when the graph partitioning is done.
01004  *
01005  * This function also returns a postorderd separator tree (CParent), and a
01006  * mapping of nodes in the graph to nodes in the separator tree (Cmember).
01007  *
01008  * workspace: Flag (nrow), Head (nrow+1), Iwork (4*nrow + (ncol if unsymmetric))
01009  *  Allocates a temporary matrix B=A*A' or B=A,
01010  *  and O(nnz(A)) temporary memory space.
01011  *  Allocates an additional 3*n*sizeof(Int) temporary workspace
01012  */
01013 
01014 UF_long CHOLMOD(nested_dissection) /* returns # of components, or -1 if error */
01015 (
01016     /* ---- input ---- */
01017     cholmod_sparse *A,  /* matrix to order */
01018     Int *fset,    /* subset of 0:(A->ncol)-1 */
01019     size_t fsize, /* size of fset */
01020     /* ---- output --- */
01021     Int *Perm,    /* size A->nrow, output permutation */
01022     Int *CParent, /* size A->nrow.  On output, CParent [c] is the parent
01023        * of component c, or EMPTY if c is a root, and where
01024        * c is in the range 0 to # of components minus 1 */
01025     Int *Cmember, /* size A->nrow.  Cmember [j] = c if node j of A is
01026        * in component c */
01027     /* --------------- */
01028     cholmod_common *Common
01029 )
01030 {
01031     double prune_dense, nd_oksep ;
01032     Int *Bp, *Bi, *Bnz, *Cstack, *Imap, *Map, *Flag, *Head, *Next, *Bnw, *Iwork,
01033   *Ipost, *NewParent, *Hash, *Cmap, *Cp, *Ci, *Cew, *Cnw, *Part, *Post,
01034   *Work3n ;
01035     unsigned Int hash ;
01036     Int n, bnz, top, i, j, k, cnode, cdense, p, cj, cn, ci, cnz, mark, c, uncol,
01037   sepsize, parent, ncomponents, threshold, ndense, pstart, pdest, pend,
01038   nd_compress, nd_camd, csize, jnext, nd_small, total_weight,
01039   nchild, child = EMPTY ;
01040     cholmod_sparse *B, *C ;
01041     size_t s ;
01042     int ok = TRUE ;
01043     DEBUG (Int cnt) ;
01044 
01045     /* ---------------------------------------------------------------------- */
01046     /* get inputs */
01047     /* ---------------------------------------------------------------------- */
01048 
01049     RETURN_IF_NULL_COMMON (EMPTY) ;
01050     RETURN_IF_NULL (A, EMPTY) ;
01051     RETURN_IF_NULL (Perm, EMPTY) ;
01052     RETURN_IF_NULL (CParent, EMPTY) ;
01053     RETURN_IF_NULL (Cmember, EMPTY) ;
01054     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
01055     Common->status = CHOLMOD_OK ;
01056 
01057     /* ---------------------------------------------------------------------- */
01058     /* quick return */
01059     /* ---------------------------------------------------------------------- */
01060 
01061     n = A->nrow ;
01062     if (n == 0)
01063     {
01064   return (1) ;
01065     }
01066 
01067     /* ---------------------------------------------------------------------- */
01068     /* get inputs */
01069     /* ---------------------------------------------------------------------- */
01070 
01071     ASSERT (CHOLMOD(dump_sparse) (A, "A for NESDIS:", Common) >= 0) ;
01072 
01073     /* get ordering parameters */
01074     prune_dense = Common->method [Common->current].prune_dense ;
01075     nd_compress = Common->method [Common->current].nd_compress ;
01076     nd_oksep = Common->method [Common->current].nd_oksep ;
01077     nd_oksep = MAX (0, nd_oksep) ;
01078     nd_oksep = MIN (1, nd_oksep) ;
01079     nd_camd = Common->method [Common->current].nd_camd ;
01080     nd_small = Common->method [Common->current].nd_small ;
01081     nd_small = MAX (4, nd_small) ;
01082 
01083     PRINT0 (("nd_components %d nd_small %d nd_oksep %g\n", 
01084   Common->method [Common->current].nd_components,
01085   nd_small, nd_oksep)) ;
01086 
01087     /* ---------------------------------------------------------------------- */
01088     /* allocate workspace */
01089     /* ---------------------------------------------------------------------- */
01090 
01091     /* s = 4*n + uncol */
01092     uncol = (A->stype == 0) ? A->ncol : 0 ;
01093     s = CHOLMOD(mult_size_t) (n, 4, &ok) ;
01094     s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
01095     if (!ok)
01096     {
01097   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
01098   return (EMPTY) ;
01099     }
01100 
01101     CHOLMOD(allocate_work) (n, s, 0, Common) ;
01102     if (Common->status < CHOLMOD_OK)
01103     {
01104   return (EMPTY) ;
01105     }
01106     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
01107 
01108     /* ---------------------------------------------------------------------- */
01109     /* get workspace */
01110     /* ---------------------------------------------------------------------- */
01111 
01112     Flag = Common->Flag ; /* size n */
01113     Head = Common->Head ; /* size n+1, all equal to -1 */
01114 
01115     Iwork = Common->Iwork ;
01116     Imap = Iwork ;    /* size n, same as Queue in find_components */
01117     Map  = Iwork + n ;    /* size n */
01118     Bnz  = Iwork + 2*((size_t) n) ; /* size n */
01119     Hash = Iwork + 3*((size_t) n) ; /* size n */
01120 
01121     Work3n = CHOLMOD(malloc) (n, 3*sizeof (Int), Common) ;
01122     Part = Work3n ;   /* size n */
01123     Bnw  = Part + n ;   /* size n */
01124     Cnw  = Bnw + n ;    /* size n */
01125 
01126     Cstack = Perm ;   /* size n, use Perm as workspace for Cstack [ */
01127     Cmap = Cmember ;    /* size n, use Cmember as workspace [ */
01128 
01129     if (Common->status < CHOLMOD_OK)
01130     {
01131   return (EMPTY) ;
01132     }
01133 
01134     /* ---------------------------------------------------------------------- */
01135     /* convert B to symmetric form with both upper/lower parts present */
01136     /* ---------------------------------------------------------------------- */
01137 
01138     /* B = A+A', A*A', or A(:,f)*A(:,f)', upper and lower parts present */
01139 
01140     if (A->stype)
01141     {
01142   /* Add the upper/lower part to a symmetric lower/upper matrix by
01143    * converting to unsymmetric mode */
01144   /* workspace: Iwork (nrow) */
01145   B = CHOLMOD(copy) (A, 0, -1, Common) ;
01146     }
01147     else
01148     {
01149   /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
01150   /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
01151   B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
01152     }
01153 
01154     if (Common->status < CHOLMOD_OK)
01155     {
01156   CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
01157   return (EMPTY) ;
01158     }
01159     Bp = B->p ;
01160     Bi = B->i ;
01161     bnz = CHOLMOD(nnz) (B, Common) ;
01162     ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;
01163     csize = MAX (n, bnz) ;
01164     ASSERT (CHOLMOD(dump_sparse) (B, "B for nd:", Common) >= 0) ;
01165 
01166     /* ---------------------------------------------------------------------- */
01167     /* initializations */
01168     /* ---------------------------------------------------------------------- */
01169 
01170     /* all nodes start out unmarked and unordered (Type 4, see below) */
01171     Common->mark = EMPTY ;
01172     CHOLMOD(clear_flag) (Common) ;
01173     ASSERT (Flag == Common->Flag) ;
01174     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
01175 
01176     for (j = 0 ; j < n ; j++)
01177     {
01178   CParent [j] = -2 ;
01179     }
01180 
01181     /* prune dense nodes from B */
01182     if (IS_NAN (prune_dense) || prune_dense < 0)
01183     {
01184   /* only remove completely dense nodes */
01185   threshold = n-2 ;
01186     }
01187     else
01188     {
01189   /* remove nodes with degree more than threshold */
01190   threshold = (Int) (MAX (16, prune_dense * sqrt ((double) (n)))) ;
01191   threshold = MIN (n, threshold) ;
01192     }
01193     ndense = 0 ;
01194     cnode = EMPTY ;
01195     cdense = EMPTY ;
01196 
01197     for (j = 0 ; j < n ; j++)
01198     {
01199   Bnz [j] = Bp [j+1] - Bp [j] ;
01200   if (Bnz [j] > threshold)
01201   {
01202       /* node j is dense, prune it from B */
01203       PRINT2 (("j is dense %d\n", j)) ;
01204       ndense++ ;
01205       if (cnode == EMPTY)
01206       {
01207     /* first dense node found becomes root of this component,
01208      * which contains all of the dense nodes found here */
01209     cdense = j ;
01210     cnode = j ;
01211     CParent [cnode] = EMPTY ;
01212       }
01213       Flag [j] = FLIP (cnode) ;
01214   }
01215     }
01216     B->packed = FALSE ;
01217     ASSERT (B->nz == NULL) ;
01218 
01219     if (ndense == n)
01220     {
01221   /* all nodes removed: Perm is identity, all nodes in component zero,
01222    * and the separator tree has just one node. */
01223   PRINT2 (("all nodes are dense\n")) ;
01224   for (k = 0 ; k < n ; k++)
01225   {
01226       Perm [k] = k ;
01227       Cmember [k] = 0 ;
01228   }
01229   CParent [0] = EMPTY ;
01230   CHOLMOD(free_sparse) (&B, Common) ;
01231   CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
01232   Common->mark = EMPTY ;
01233   CHOLMOD(clear_flag) (Common) ;
01234   return (1) ;
01235     }
01236 
01237     /* Cp and Ci are workspace to construct the subgraphs to partition */
01238     C = CHOLMOD(allocate_sparse) (n, n, csize, FALSE, TRUE, 0, CHOLMOD_PATTERN,
01239       Common) ;
01240     Cew  = CHOLMOD(malloc) (csize, sizeof (Int), Common) ;
01241 
01242     if (Common->status < CHOLMOD_OK)
01243     {
01244   /* out of memory */
01245   CHOLMOD(free_sparse) (&C, Common) ;
01246   CHOLMOD(free_sparse) (&B, Common) ;
01247   CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
01248   CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
01249   Common->mark = EMPTY ;
01250   CHOLMOD(clear_flag) (Common) ;
01251   PRINT2 (("out of memory for C, etc\n")) ;
01252   return (EMPTY) ;
01253     }
01254 
01255     Cp = C->p ;
01256     Ci = C->i ;
01257 
01258     /* create initial unit node and edge weights */
01259     for (j = 0 ; j < n ; j++)
01260     {
01261   Bnw [j] = 1 ;
01262     }
01263     for (p = 0 ; p < csize ; p++)
01264     {
01265   Cew [p] = 1 ;
01266     }
01267 
01268     /* push the initial connnected components of B onto the Cstack */
01269     top = EMPTY ; /* Cstack is empty */
01270     /* workspace: Flag (nrow), Iwork (nrow); use Imap as workspace for Queue [*/
01271     find_components (B, NULL, n, cnode, NULL,
01272       Bnz, CParent, Cstack, &top, Imap, Common) ;
01273     /* done using Imap as workspace for Queue ] */
01274 
01275     /* Nodes can now be of Type 0, 1, 2, or 4 (see definition below) */
01276 
01277     /* ---------------------------------------------------------------------- */
01278     /* while Cstack is not empty, do: */
01279     /* ---------------------------------------------------------------------- */
01280 
01281     while (top >= 0)
01282     {
01283 
01284   /* clear the Flag array, but do not modify negative entries in Flag  */
01285   mark = clear_flag (Common) ;
01286 
01287   DEBUG (for (i = 0 ; i < n ; i++) Imap [i] = EMPTY) ;
01288 
01289   /* ------------------------------------------------------------------ */
01290   /* get node(s) from the top of the Cstack */
01291   /* ------------------------------------------------------------------ */
01292 
01293   /* i is the repnode of its (unordered) connected component.  Get
01294    * all repnodes for all connected components of a single part.  If
01295    * each connected component is to be ordered separately (nd_components
01296    * is TRUE), then this while loop iterates just once. */
01297 
01298   cnode = EMPTY ;
01299   cn = 0 ;
01300   while (cnode == EMPTY)
01301   {
01302       i = Cstack [top--] ;
01303 
01304       if (i < 0)
01305       {
01306     /* this is the last node in this component */
01307     i = FLIP (i) ;
01308     cnode = i ;
01309       }
01310 
01311       ASSERT (i >= 0 && i < n && Flag [i] >= EMPTY) ;
01312 
01313       /* place i in the queue and mark it */
01314       Map [cn] = i ;
01315       Flag [i] = mark ;
01316       Imap [i] = cn ;
01317       cn++ ;
01318   }
01319 
01320   ASSERT (cnode != EMPTY) ;
01321 
01322   /* During ordering, there are five kinds of nodes in the graph of B,
01323    * based on Flag [j] and CParent [j] for nodes j = 0 to n-1:
01324    *
01325    * Type 0: If cnode is a repnode of an unordered component, then
01326    * CParent [cnode] is in the range EMPTY to n-1 and
01327    * Flag [cnode] >= EMPTY.  This is a "live" node.
01328    *
01329    * Type 1: If cnode is a repnode of an ordered separator component,
01330    * then Flag [cnode] < EMPTY and FLAG [cnode] = FLIP (cnode).
01331    * CParent [cnode] is in the range EMPTY to n-1.  cnode is a root of
01332    * the separator tree if CParent [cnode] == EMPTY.  This node is dead.
01333    *
01334    * Type 2: If node j isn't a repnode, has not been absorbed via
01335    * graph compression into another node, but is in an ordered separator
01336    * component, then cnode = FLIP (Flag [j]) gives the repnode of the
01337    * component that contains j and CParent [j]  is -2.  This node is dead.
01338    * Note that Flag [j] < EMPTY.
01339    *
01340    * Type 3: If node i has been absorbed via graph compression into some
01341    * other node j = FLIP (Flag [i]) where j is not a repnode.
01342    * CParent [j] is -2.  Node i may or may not be in an ordered
01343    * component.  This node is dead.  Note that Flag [j] < EMPTY.
01344    *
01345    * Type 4: If node j is "live" (not in an ordered component, and not
01346    * absorbed into any other node), then Flag [j] >= EMPTY.
01347    *
01348    * Only "live" nodes (of type 0 or 4) are placed in a subgraph to be
01349    * partitioned.  Node j is alive if Flag [j] >= EMPTY, and dead if
01350    * Flag [j] < EMPTY.
01351    */
01352 
01353   /* ------------------------------------------------------------------ */
01354   /* create the subgraph for this connected component C */
01355   /* ------------------------------------------------------------------ */
01356 
01357   /* Do a breadth-first search of the graph starting at cnode.
01358    * use Map [0..cn-1] for nodes in the component C [
01359    * use Cnw and Cew for node and edge weights of the resulting subgraph [
01360    * use Cp and Ci for the resulting subgraph [
01361    * use Imap [i] for all nodes i in B that are in the component C [
01362    */
01363 
01364   cnz = 0 ;
01365   total_weight = 0 ;
01366   for (cj = 0 ; cj < cn ; cj++)
01367   {
01368       /* get node j from the head of the queue; it is node cj of C */
01369       j = Map [cj] ;
01370       ASSERT (Flag [j] == mark) ;
01371       Cp [cj] = cnz ;
01372       Cnw [cj] = Bnw [j] ;
01373       ASSERT (Cnw [cj] >= 0) ;
01374       total_weight += Cnw [cj] ;
01375       pstart = Bp [j] ;
01376       pdest = pstart ;
01377       pend = pstart + Bnz [j] ;
01378       hash = cj ;
01379       for (p = pstart ; p < pend ; p++)
01380       {
01381     i = Bi [p] ;
01382     /* prune diagonal entries and dead edges from B */
01383     if (i != j && Flag [i] >= EMPTY)
01384     {
01385         /* live node i is in the current component */
01386         Bi [pdest++] = i ;
01387         if (Flag [i] != mark)
01388         {
01389       /* First time node i has been seen, it is a new node
01390        * of C.  place node i in the queue and mark it */
01391       Map [cn] = i ;
01392       Flag [i] = mark ;
01393       Imap [i] = cn ;
01394       cn++ ;
01395         }
01396         /* place the edge (cj,ci) in the adjacency list of cj */
01397         ci = Imap [i] ;
01398         ASSERT (ci >= 0 && ci < cn && ci != cj && cnz < csize) ;
01399         Ci [cnz++] = ci ;
01400         hash += ci ;
01401     }
01402       }
01403       /* edges to dead nodes have been removed */
01404       Bnz [j] = pdest - pstart ;
01405       /* finalize the hash key for column j */
01406       hash %= csize ;
01407       Hash [cj] = (Int) hash ;
01408       ASSERT (Hash [cj] >= 0 && Hash [cj] < csize) ;
01409   }
01410   Cp [cn] = cnz ;
01411   C->nrow = cn ;
01412   C->ncol = cn ;  /* affects mem stats unless restored when C free'd */
01413 
01414   /* contents of Imap no longer needed ] */
01415 
01416 #ifndef NDEBUG
01417   for (cj = 0 ; cj < cn ; cj++)
01418   {
01419       j = Map [cj] ;
01420       PRINT2 (("----------------------------C column cj: "ID" j: "ID"\n",
01421     cj, j)) ;
01422       ASSERT (j >= 0 && j < n) ;
01423       ASSERT (Flag [j] >= EMPTY) ;
01424       for (p = Cp [cj] ; p < Cp [cj+1] ; p++)
01425       {
01426     ci = Ci [p] ;
01427     i = Map [ci] ;
01428     PRINT3 (("ci: "ID" i: "ID"\n", ci, i)) ;
01429     ASSERT (ci != cj && ci >= 0 && ci < cn) ;
01430     ASSERT (i != j && i >= 0 && i < n) ;
01431     ASSERT (Flag [i] >= EMPTY) ;
01432       }
01433   }
01434 #endif
01435 
01436   PRINT0 (("consider cn %d nd_small %d ", cn, nd_small)) ;
01437   if (cn < nd_small)  /* could be 'total_weight < nd_small' instead */
01438   {
01439       /* place all nodes in the separator */
01440       PRINT0 ((" too small\n")) ;
01441       sepsize = total_weight ;
01442   }
01443   else
01444   {
01445 
01446       /* Cp and Ci now contain the component, with cn nodes and cnz
01447        * nonzeros.  The mapping of a node cj into node j the main graph
01448        * B is given by Map [cj] = j */
01449       PRINT0 ((" cut\n")) ;
01450 
01451       /* -------------------------------------------------------------- */
01452       /* compress and partition the graph C */
01453       /* -------------------------------------------------------------- */
01454 
01455       /* The edge weights Cew [0..csize-1] are all 1's on input to and
01456        * output from the partition routine. */
01457 
01458       sepsize = partition (
01459 #ifndef NDEBUG
01460         csize,
01461 #endif
01462         nd_compress, Hash, C, Cnw, Cew,
01463         Cmap, Part, Common) ;
01464 
01465       /* contents of Cp and Ci no longer needed ] */
01466 
01467       if (sepsize < 0)
01468       {
01469     /* failed */
01470     C->ncol = n ;   /* restore size for memory usage statistics */
01471     CHOLMOD(free_sparse) (&C, Common) ;
01472     CHOLMOD(free_sparse) (&B, Common) ;
01473     CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
01474     CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
01475     Common->mark = EMPTY ;
01476     CHOLMOD(clear_flag) (Common) ;
01477     return (EMPTY) ;
01478       }
01479 
01480       /* -------------------------------------------------------------- */
01481       /* compress B based on how C was compressed */
01482       /* -------------------------------------------------------------- */
01483 
01484       for (ci = 0 ; ci < cn ; ci++)
01485       {
01486     if (Hash [ci] < EMPTY)
01487     {
01488         /* ci is dead in C, having been absorbed into cj */
01489         cj = FLIP (Hash [ci]) ;
01490         PRINT2 (("In C, "ID" absorbed into "ID" (wgt now "ID")\n",
01491           ci, cj, Cnw [cj])) ;
01492         /* i is dead in B, having been absorbed into j */
01493         i = Map [ci] ;
01494         j = Map [cj] ;
01495         PRINT2 (("In B, "ID" (wgt "ID") => "ID" (wgt "ID")\n",
01496         i, Bnw [i], j, Bnw [j], Cnw [cj])) ;
01497         /* more than one node may be absorbed into j.  This is
01498          * accounted for in Cnw [cj].  Assign it here rather
01499          * than += Bnw [i] */
01500         Bnw [i] = 0 ;
01501         Bnw [j] = Cnw [cj] ;
01502         Flag [i] = FLIP (j) ;
01503     }
01504       }
01505 
01506       DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) cnt += Bnw [j]) ;
01507       ASSERT (cnt == n) ;
01508   }
01509 
01510   /* contents of Cnw [0..cn-1] no longer needed ] */
01511 
01512   /* ------------------------------------------------------------------ */
01513   /* order the separator, and stack the components when C is split */
01514   /* ------------------------------------------------------------------ */
01515 
01516   /* one more component has been found: either the separator of C,
01517    * or all of C */
01518 
01519   ASSERT (sepsize >= 0 && sepsize <= total_weight) ;
01520 
01521   PRINT0 (("sepsize %d tot %d : %8.4f ", sepsize, total_weight,
01522       ((double) sepsize) / ((double) total_weight))) ;
01523 
01524   if (sepsize == total_weight || sepsize == 0 ||
01525       sepsize > nd_oksep * total_weight)
01526   {
01527       /* Order the nodes in the component.  The separator is too large,
01528        * or empty.  Note that the partition routine cannot return a
01529        * sepsize of zero, but it can return a separator consisting of the
01530        * whole graph.  The "sepsize == 0" test is kept, above, in case the
01531        * partition routine changes.  In either case, this component
01532        * remains unsplit, and becomes a leaf of the separator tree. */
01533       PRINT2 (("cnode %d sepsize zero or all of graph: "ID"\n",
01534     cnode, sepsize)) ;
01535       for (cj = 0 ; cj < cn ; cj++)
01536       {
01537     j = Map [cj] ;
01538     Flag [j] = FLIP (cnode) ;
01539     PRINT2 (("      node cj: "ID" j: "ID" ordered\n", cj, j)) ;
01540       }
01541       ASSERT (Flag [cnode] == FLIP (cnode)) ;
01542       ASSERT (cnode != EMPTY && Flag [cnode] < EMPTY) ;
01543       PRINT0 (("discarded\n")) ;
01544 
01545   }
01546   else
01547   {
01548 
01549       /* Order the nodes in the separator of C and find a new repnode
01550        * cnode that is in the separator of C.  This requires the separator
01551        * to be non-empty. */
01552       PRINT0 (("sepsize not tiny: "ID"\n", sepsize)) ;
01553       parent = CParent [cnode] ;
01554       ASSERT (parent >= EMPTY && parent < n) ;
01555       CParent [cnode] = -2 ;
01556       cnode = EMPTY ;
01557       for (cj = 0 ; cj < cn ; cj++)
01558       {
01559     j = Map [cj] ;
01560     if (Part [cj] == 2)
01561     {
01562         /* All nodes in the separator become part of a component
01563          * whose repnode is cnode */
01564         PRINT2 (("node cj: "ID" j: "ID" ordered\n", cj, j)) ;
01565         if (cnode == EMPTY)
01566         {
01567       PRINT2(("------------new cnode: cj "ID" j "ID"\n",
01568             cj, j)) ;
01569       cnode = j ;
01570         }
01571         Flag [j] = FLIP (cnode) ;
01572     }
01573     else
01574     {
01575         PRINT2 (("      node cj: "ID" j: "ID" not ordered\n",
01576         cj, j)) ;
01577     }
01578       }
01579       ASSERT (cnode != EMPTY && Flag [cnode] < EMPTY) ;
01580       ASSERT (CParent [cnode] == -2) ;
01581       CParent [cnode] = parent ;
01582 
01583       /* find the connected components when C is split, and push
01584        * them on the Cstack.  Use Imap as workspace for Queue. [ */
01585       /* workspace: Flag (nrow) */
01586       find_components (B, Map, cn, cnode, Part, Bnz,
01587         CParent, Cstack, &top, Imap, Common) ;
01588       /* done using Imap as workspace for Queue ] */
01589   }
01590   /* contents of Map [0..cn-1] no longer needed ] */
01591     }
01592 
01593     /* done using Cmember as workspace for Cmap ] */
01594     /* done using Perm as workspace for Cstack ] */
01595 
01596     /* ---------------------------------------------------------------------- */
01597     /* place nodes removed via compression into their proper component */
01598     /* ---------------------------------------------------------------------- */
01599 
01600     /* At this point, all nodes are of Type 1, 2, or 3, as defined above. */
01601 
01602     for (i = 0 ; i < n ; i++)
01603     {
01604   /* find the repnode cnode that contains node i */
01605   j = FLIP (Flag [i]) ;
01606   PRINT2 (("\nfind component for "ID", in: "ID"\n", i, j)) ;
01607   ASSERT (j >= 0 && j < n) ;
01608   DEBUG (cnt = 0) ;
01609   while (CParent [j] == -2)
01610   {
01611       j = FLIP (Flag [j]) ;
01612       PRINT2 (("    walk up to "ID" ", j)) ;
01613       ASSERT (j >= 0 && j < n) ;
01614       PRINT2 ((" CParent "ID"\n", CParent [j])) ;
01615       ASSERT (cnt < n) ;
01616       DEBUG (cnt++) ;
01617   }
01618   cnode = j ;
01619   ASSERT (cnode >= 0 && cnode < n) ;
01620   ASSERT (CParent [cnode] >= EMPTY && CParent [cnode] < n) ;
01621   PRINT2 (("i "ID" is in component with cnode "ID"\n", i, cnode)) ;
01622   ASSERT (Flag [cnode] == FLIP (cnode)) ;
01623 
01624   /* Mark all nodes along the path from i to cnode as being in the
01625    * component whos repnode is cnode.  Perform path compression.  */
01626   j = FLIP (Flag [i]) ;
01627   Flag [i] = FLIP (cnode) ;
01628   DEBUG (cnt = 0) ;
01629   while (CParent [j] == -2)
01630   {
01631       ASSERT (j >= 0 && j < n) ;
01632       jnext = FLIP (Flag [j]) ;
01633       PRINT2 (("    "ID" walk "ID" set cnode to "ID"\n", i, j, cnode)) ;
01634       ASSERT (cnt < n) ;
01635       DEBUG (cnt++) ;
01636       Flag [j] = FLIP (cnode) ;
01637       j = jnext ;
01638   }
01639     }
01640 
01641     /* At this point, all nodes fall into Types 1 or 2, as defined above. */
01642 
01643 #ifndef NDEBUG
01644     for (j = 0 ; j < n ; j++)
01645     {
01646   PRINT2 (("j %d CParent %d  ", j, CParent [j])) ;
01647   if (CParent [j] >= EMPTY && CParent [j] < n)
01648   {
01649       /* case 1: j is a repnode of a component */
01650       cnode = j ;
01651       PRINT2 ((" a repnode\n")) ;
01652   }
01653   else
01654   {
01655       /* case 2: j is not a repnode of a component */
01656       cnode = FLIP (Flag [j]) ;
01657       PRINT2 ((" repnode is %d\n", cnode)) ;
01658       ASSERT (cnode >= 0 && cnode < n) ;
01659       ASSERT (CParent [cnode] >= EMPTY && CParent [cnode] < n) ;
01660   }
01661   ASSERT (Flag [cnode] == FLIP (cnode)) ;
01662   /* case 3 no longer holds */
01663     }
01664 #endif
01665 
01666     /* ---------------------------------------------------------------------- */
01667     /* free workspace */
01668     /* ---------------------------------------------------------------------- */
01669 
01670     C->ncol = n ;   /* restore size for memory usage statistics */
01671     CHOLMOD(free_sparse) (&C, Common) ;
01672     CHOLMOD(free_sparse) (&B, Common) ;
01673     CHOLMOD(free) (csize, sizeof (Int), Cew, Common) ;
01674     CHOLMOD(free) (3*n, sizeof (Int), Work3n, Common) ;
01675 
01676     /* ---------------------------------------------------------------------- */
01677     /* handle dense nodes */
01678     /* ---------------------------------------------------------------------- */
01679 
01680     /* The separator tree has nodes with either no children or two or more
01681      * children - with one exception.  There may exist a single root node with
01682      * exactly one child, which holds the dense rows/columns of the matrix.
01683      * Delete this node if it exists. */
01684 
01685     if (ndense > 0)
01686     {
01687   ASSERT (CParent [cdense] == EMPTY) ;  /* cdense has no parent */
01688   /* find the children of cdense */
01689   nchild = 0 ;
01690   for (j = 0 ; j < n ; j++)
01691   {
01692       if (CParent [j] == cdense)
01693       {
01694     nchild++ ;
01695     child = j ;
01696       }
01697   }
01698   if (nchild == 1)
01699   {
01700       /* the cdense node has just one child; merge the two nodes */
01701       PRINT1 (("root has one child\n")) ;
01702       CParent [cdense] = -2 ;   /* cdense is deleted */
01703       CParent [child] = EMPTY ;   /* child becomes a root */
01704       for (j = 0 ; j < n ; j++)
01705       {
01706     if (Flag [j] == FLIP (cdense))
01707     {
01708         /* j is a dense node */
01709         PRINT1 (("dense %d\n", j)) ;
01710         Flag [j] = FLIP (child) ;
01711     }
01712       }
01713   }
01714     }
01715 
01716     /* ---------------------------------------------------------------------- */
01717     /* postorder the components */
01718     /* ---------------------------------------------------------------------- */
01719 
01720     DEBUG (for (cnt = 0, j = 0 ; j < n ; j++) if (CParent [j] != -2) cnt++) ;
01721 
01722     /* use Cmember as workspace for Post [ */
01723     Post = Cmember ;
01724 
01725     /* cholmod_postorder uses Head and Iwork [0..2n].  It does not use Flag,
01726      * which here holds the mapping of nodes to repnodes.  It ignores all nodes
01727      * for which CParent [j] < -1, so it operates just on the repnodes. */
01728     /* workspace: Head (n), Iwork (2*n) */
01729     ncomponents = CHOLMOD(postorder) (CParent, n, NULL, Post, Common) ;
01730     ASSERT (cnt == ncomponents) ;
01731 
01732     /* use Iwork [0..n-1] as workspace for Ipost ( */
01733     Ipost = Iwork ;
01734     DEBUG (for (j = 0 ; j < n ; j++) Ipost [j] = EMPTY) ;
01735 
01736     /* compute inverse postorder */
01737     for (c = 0 ; c < ncomponents ; c++)
01738     {
01739   cnode = Post [c] ;
01740   ASSERT (cnode >= 0 && cnode < n) ;
01741   Ipost [cnode] = c ;
01742   ASSERT (Head [c] == EMPTY) ;
01743     }
01744 
01745     /* adjust the parent array */
01746     /* Iwork [n..2n-1] used for NewParent [ */
01747     NewParent = Iwork + n ;
01748     for (c = 0 ; c < ncomponents ; c++)
01749     {
01750   parent = CParent [Post [c]] ;
01751   NewParent [c] = (parent == EMPTY) ? EMPTY : (Ipost [parent]) ;
01752     }
01753     for (c = 0 ; c < ncomponents ; c++)
01754     {
01755   CParent [c] = NewParent [c] ;
01756     }
01757     ASSERT (CHOLMOD(dump_parent) (CParent, ncomponents, "CParent", Common)) ;
01758 
01759     /* Iwork [n..2n-1] no longer needed for NewParent ] */
01760     /* Cmember no longer needed for Post ] */
01761 
01762 #ifndef NDEBUG
01763     /* count the number of children of each node */
01764     for (c = 0 ; c < ncomponents ; c++)
01765     {
01766   Cmember [c] = 0 ;
01767     }
01768     for (c = 0 ; c < ncomponents ; c++)
01769     {
01770   if (CParent [c] != EMPTY) Cmember [CParent [c]]++ ;
01771     }
01772     for (c = 0 ; c < ncomponents ; c++)
01773     {
01774   /* a node is either a leaf, or has 2 or more children */
01775   ASSERT (Cmember [c] == 0 || Cmember [c] >= 2) ;
01776     }
01777 #endif
01778 
01779     /* ---------------------------------------------------------------------- */
01780     /* place each node in its component */
01781     /* ---------------------------------------------------------------------- */
01782 
01783     for (j = 0 ; j < n ; j++)
01784     {
01785   /* node j is in the cth component, whose repnode is cnode */
01786   cnode = FLIP (Flag [j]) ;
01787   PRINT2 (("j "ID"  flag "ID" cnode "ID"\n",
01788         j, Flag [j], FLIP (Flag [j]))) ;
01789   ASSERT (cnode >= 0 && cnode < n) ;
01790   c = Ipost [cnode] ;
01791   ASSERT (c >= 0 && c < ncomponents) ;
01792   Cmember [j] = c ;
01793     }
01794 
01795     /* Flag no longer needed for the node-to-component mapping */
01796 
01797     /* done using Iwork [0..n-1] as workspace for Ipost ) */
01798 
01799     /* ---------------------------------------------------------------------- */
01800     /* clear the Flag array */
01801     /* ---------------------------------------------------------------------- */
01802 
01803     Common->mark = EMPTY ;
01804     CHOLMOD(clear_flag) (Common) ;
01805     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
01806 
01807     /* ---------------------------------------------------------------------- */
01808     /* find the permutation */
01809     /* ---------------------------------------------------------------------- */
01810 
01811     PRINT1 (("nd_camd: %d A->stype %d\n", nd_camd, A->stype)) ;
01812 
01813     if (nd_camd)
01814     {
01815 
01816   /* ------------------------------------------------------------------ */
01817   /* apply camd, csymamd, or ccolamd using the Cmember constraints */
01818   /* ------------------------------------------------------------------ */
01819 
01820   if (A->stype != 0)
01821   {
01822       /* ordering A+A', so fset and fsize are ignored.
01823        * Add the upper/lower part to a symmetric lower/upper matrix by
01824        * converting to unsymmetric mode
01825        * workspace: Iwork (nrow) */
01826       B = CHOLMOD(copy) (A, 0, -1, Common) ;
01827       if (Common->status < CHOLMOD_OK)
01828       {
01829     PRINT0 (("make symmetric failed\n")) ;
01830     return (EMPTY) ;
01831       }
01832       ASSERT ((Int) (B->nrow) == n && (Int) (B->ncol) == n) ;
01833       PRINT2 (("nested dissection (2)\n")) ;
01834       B->stype = -1 ;
01835       if (nd_camd == 2)
01836       {
01837     /* workspace:  Head (nrow+1), Iwork (nrow) if symmetric-upper */
01838     ok = CHOLMOD(csymamd) (B, Cmember, Perm, Common) ;
01839       }
01840       else
01841       {
01842     /* workspace: Head (nrow), Iwork (4*nrow) */
01843     ok = CHOLMOD(camd) (B, NULL, 0, Cmember, Perm, Common) ;
01844       }
01845       CHOLMOD(free_sparse) (&B, Common) ;
01846       if (!ok)
01847       {
01848     /* failed */
01849     PRINT0 (("camd/csymamd failed\n")) ;
01850     return (EMPTY) ;
01851       }
01852   }
01853   else
01854   {
01855       /* ordering A*A' or A(:,f)*A(:,f)' */
01856       /* workspace: Iwork (nrow if no fset; MAX(nrow,ncol) if fset) */
01857       if (!CHOLMOD(ccolamd) (A, fset, fsize, Cmember, Perm, Common))
01858       {
01859     /* ccolamd failed */
01860     PRINT2 (("ccolamd failed\n")) ;
01861     return (EMPTY) ;
01862       }
01863   }
01864 
01865     }
01866     else
01867     {
01868 
01869   /* ------------------------------------------------------------------ */
01870   /* natural ordering of each component */
01871   /* ------------------------------------------------------------------ */
01872 
01873   /* use Iwork [0..n-1] for Next [ */
01874   Next = Iwork  ;
01875 
01876   /* ------------------------------------------------------------------ */
01877   /* place the nodes in link lists, one list per component */
01878   /* ------------------------------------------------------------------ */
01879 
01880   /* do so in reverse order, to preserve original ordering */
01881   for (j = n-1 ; j >= 0 ; j--)
01882   {
01883       /* node j is in the cth component */
01884       c = Cmember [j] ;
01885       ASSERT (c >= 0 && c < ncomponents) ;
01886       /* place node j in link list for component c */
01887       Next [j] = Head [c] ;
01888       Head [c] = j ;
01889   }
01890 
01891   /* ------------------------------------------------------------------ */
01892   /* order each node in each component */
01893   /* ------------------------------------------------------------------ */
01894 
01895   k = 0 ;
01896   for (c = 0 ; c < ncomponents ; c++)
01897   {
01898       for (j = Head [c] ; j != EMPTY ; j = Next [j])
01899       {
01900     Perm [k++] = j ;
01901       }
01902       Head [c] = EMPTY ;
01903   }
01904   ASSERT (k == n) ;
01905 
01906   /* done using Iwork [0..n-1] for Next ] */
01907     }
01908 
01909     /* ---------------------------------------------------------------------- */
01910     /* clear workspace and return number of components */
01911     /* ---------------------------------------------------------------------- */
01912 
01913     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
01914     return (ncomponents) ;
01915 }
01916 
01917 /* ========================================================================== */
01918 /* === cholmod_collapse_septree ============================================= */
01919 /* ========================================================================== */
01920 
01921 /* cholmod_nested_dissection returns the separator tree that was used in the
01922  * constrained minimum degree algorithm.  Parameter settings (nd_small,
01923  * nd_oksep, etc) that give a good fill-reducing ordering may give too fine of
01924  * a separator tree for other uses (parallelism, multi-level LPDASA, etc).  This
01925  * function takes as input the separator tree computed by
01926  * cholmod_nested_dissection, and collapses selected subtrees into single
01927  * nodes.  A subtree is collapsed if its root node (the separator) is large
01928  * compared to the total number of nodes in the subtree, or if the subtree is
01929  * small.  Note that the separator tree may actually be a forest.
01930  *
01931  * nd_oksep and nd_small act just like the ordering parameters in Common.
01932  * Returns the new number of nodes in the separator tree.
01933  */
01934 
01935 UF_long CHOLMOD(collapse_septree)
01936 (
01937     /* ---- input ---- */
01938     size_t n,   /* # of nodes in the graph */
01939     size_t ncomponents, /* # of nodes in the separator tree (must be <= n) */
01940     double nd_oksep,    /* collapse if #sep >= nd_oksep * #nodes in subtree */
01941     size_t nd_small,    /* collapse if #nodes in subtree < nd_small */
01942     /* ---- in/out --- */
01943     Int *CParent, /* size ncomponents; from cholmod_nested_dissection */
01944     Int *Cmember, /* size n; from cholmod_nested_dissection */
01945     /* --------------- */
01946     cholmod_common *Common
01947 )
01948 {
01949     Int *First, *Count, *Csubtree, *W, *Map ;
01950     Int c, j, k, nc, sepsize, total_weight, parent, nc_new, first ;
01951     int collapse = FALSE, ok = TRUE ;
01952     size_t s ;
01953 
01954     /* ---------------------------------------------------------------------- */
01955     /* get inputs */
01956     /* ---------------------------------------------------------------------- */
01957 
01958     RETURN_IF_NULL_COMMON (EMPTY) ;
01959     RETURN_IF_NULL (CParent, EMPTY) ;
01960     RETURN_IF_NULL (Cmember, EMPTY) ;
01961     if (n < ncomponents)
01962     {
01963   ERROR (CHOLMOD_INVALID, "invalid separator tree") ;
01964   return (EMPTY) ;
01965     }
01966     Common->status = CHOLMOD_OK ;
01967     nc = ncomponents ;
01968     if (n <= 1 || ncomponents <= 1)
01969     {
01970   /* no change; tree is one node already */
01971   return (nc) ;
01972     }
01973 
01974     nd_oksep = MAX (0, nd_oksep) ;
01975     nd_oksep = MIN (1, nd_oksep) ;
01976     nd_small = MAX (4, nd_small) ;
01977 
01978     /* ---------------------------------------------------------------------- */
01979     /* allocate workspace */
01980     /* ---------------------------------------------------------------------- */
01981 
01982     /* s = 3*ncomponents */
01983     s = CHOLMOD(mult_size_t) (ncomponents, 3, &ok) ;
01984     if (!ok)
01985     {
01986   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
01987   return (EMPTY) ;
01988     }
01989     CHOLMOD(allocate_work) (0, s, 0, Common) ;
01990     if (Common->status < CHOLMOD_OK)
01991     {
01992   return (EMPTY) ;
01993     }
01994     W = Common->Iwork ;
01995     Count    = W ; W += ncomponents ;     /* size ncomponents */
01996     Csubtree = W ; W += ncomponents ;     /* size ncomponents */
01997     First    = W ; W += ncomponents ;     /* size ncomponents */
01998 
01999     /* ---------------------------------------------------------------------- */
02000     /* find the first descendant of each node of the separator tree */
02001     /* ---------------------------------------------------------------------- */
02002 
02003     for (c = 0 ; c < nc ; c++)
02004     {
02005   First [c] = EMPTY ;
02006     }
02007     for (k = 0 ; k < nc ; k++)
02008     {
02009   for (c = k ; c != EMPTY && First [c] == -1 ; c = CParent [c])
02010   {
02011       ASSERT (c >= 0 && c < nc) ;
02012       First [c] = k ;
02013   }
02014     }
02015 
02016     /* ---------------------------------------------------------------------- */
02017     /* find the number of nodes of the graph in each node of the tree */
02018     /* ---------------------------------------------------------------------- */
02019 
02020     for (c = 0 ; c < nc ; c++)
02021     {
02022   Count [c] = 0 ;
02023     }
02024     for (j = 0 ; j < (Int) n ; j++)
02025     {
02026   ASSERT (Cmember [j] >= 0 && Cmember [j] < nc) ;
02027   Count [Cmember [j]]++ ;
02028     }
02029 
02030     /* ---------------------------------------------------------------------- */
02031     /* find the number of nodes in each subtree */
02032     /* ---------------------------------------------------------------------- */
02033 
02034     for (c = 0 ; c < nc ; c++)
02035     {
02036   /* each subtree includes its root */
02037   Csubtree [c] = Count [c] ;
02038   PRINT1 ((ID" size "ID" parent "ID" first "ID"\n",
02039       c, Count [c], CParent [c], First [c])) ;
02040     }
02041 
02042     for (c = 0 ; c < nc ; c++)
02043     {
02044   /* add the subtree of the child, c, into the count of its parent */
02045   parent = CParent [c] ;
02046   ASSERT (parent >= EMPTY && parent < nc) ;
02047   if (parent != EMPTY)
02048   {
02049       Csubtree [parent] += Csubtree [c] ;
02050   }
02051     }
02052 
02053 #ifndef NDEBUG
02054     /* the sum of the roots should be n */
02055     j = 0 ;
02056     for (c = 0 ; c < nc ; c++) if (CParent [c] == EMPTY) j += Csubtree [c] ;
02057     ASSERT (j == (Int) n) ;
02058 #endif
02059 
02060     /* ---------------------------------------------------------------------- */
02061     /* find subtrees to collapse */
02062     /* ---------------------------------------------------------------------- */
02063 
02064     /* consider all nodes in reverse post-order */
02065     for (c = nc-1 ; c >= 0 ; c--)
02066     {
02067   /* consider the subtree rooted at node c */
02068   sepsize = Count [c] ;
02069   total_weight = Csubtree [c] ;
02070   PRINT1 (("Node "ID" sepsize "ID" subtree "ID" ratio %g\n", c, sepsize,
02071       total_weight, ((double) sepsize)/((double) total_weight))) ;
02072   first = First [c] ;
02073   if (first < c &&    /* c must not be a leaf */
02074      (sepsize > nd_oksep * total_weight || total_weight < (int) nd_small))
02075   {
02076       /* this separator is too large, or the subtree is too small.
02077        * collapse the tree, by converting the entire subtree rooted at
02078        * c into a single node.  The subtree consists of all nodes from
02079        * First[c] to the root c.  Flag all nodes from First[c] to c-1
02080        * as dead.
02081        */
02082       collapse = TRUE ;
02083       for (k = first ; k < c ; k++)
02084       {
02085     CParent [k] = -2 ;
02086     PRINT1 (("   collapse node "ID"\n", k)) ;
02087       }
02088       /* continue at the next node, first-1 */
02089       c = first ;
02090   }
02091     }
02092 
02093     PRINT1 (("collapse: %d\n", collapse)) ;
02094 
02095     /* ---------------------------------------------------------------------- */
02096     /* compress the tree */
02097     /* ---------------------------------------------------------------------- */
02098 
02099     Map = Count ; /* Count no longer needed */
02100 
02101     nc_new = nc ;
02102     if (collapse)
02103     {
02104   nc_new = 0 ;
02105   for (c = 0 ; c < nc ; c++)
02106   {
02107       Map [c] = nc_new ;
02108       if (CParent [c] >= EMPTY)
02109       {
02110     /* node c is alive, and becomes node Map[c] in the new tree.
02111      * Increment nc_new for the next node c. */
02112     nc_new++ ;
02113       }
02114   }
02115   PRINT1 (("Collapse the tree from "ID" to "ID" nodes\n", nc, nc_new)) ;
02116   ASSERT (nc_new > 0) ;
02117   for (c = 0 ; c < nc ; c++)
02118   {
02119       parent = CParent [c] ;
02120       if (parent >= EMPTY)
02121       {
02122     /* node c is alive */
02123     CParent [Map [c]] = (parent == EMPTY) ? EMPTY : Map [parent] ;
02124       }
02125   }
02126   for (j = 0 ; j < (Int) n ; j++)
02127   {
02128       PRINT1 (("j "ID" Cmember[j] "ID" Map[Cmember[j]] "ID"\n",
02129     j, Cmember [j], Map [Cmember [j]])) ;
02130       Cmember [j] = Map [Cmember [j]] ;
02131   }
02132     }
02133 
02134     /* ---------------------------------------------------------------------- */
02135     /* return new size of separator tree */
02136     /* ---------------------------------------------------------------------- */
02137 
02138     return (nc_new) ;
02139 }
02140 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines