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