Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_metis.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Partition/cholmod_metis ============================================== */
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 interface to the METIS package (Version 4.0.1):
00015  *
00016  * cholmod_metis_bisector:
00017  *
00018  *  Wrapper for METIS_NodeComputeSeparator.  Finds a set of nodes that
00019  *  partitions the graph into two parts.
00020  *
00021  * cholmod_metis:
00022  *
00023  *  Wrapper for METIS_NodeND, METIS's own nested dissection algorithm.
00024  *  Typically faster than cholmod_nested_dissection, mostly because it
00025  *  uses minimum degree on just the leaves of the separator tree, rather
00026  *  than the whole matrix.
00027  *
00028  * Note that METIS does not return an error if it runs out of memory.  Instead,
00029  * it terminates the program.  This interface attempts to avoid that problem
00030  * by preallocating space that should be large enough for any memory allocations
00031  * within METIS, and then freeing that space, just before the call to METIS.
00032  * While this is not guaranteed to work, it is very unlikely to fail.  If you
00033  * encounter this problem, increase Common->metis_memory.  If you don't mind
00034  * having your program terminated, set Common->metis_memory to zero (a value of
00035  * 2.0 is usually safe).  Several other METIS workarounds are made in the
00036  * routines in this file.  See the description of metis_memory_ok, just below,
00037  * for more details.
00038  *
00039  * FUTURE WORK: interfaces to other partitioners (CHACO, SCOTCH, JOSTLE, ... )
00040  *
00041  * workspace: several size-nz and size-n temporary arrays.  Uses no workspace
00042  * in Common.
00043  *
00044  * Supports any xtype (pattern, real, complex, or zomplex).
00045  */
00046 
00047 #ifndef NPARTITION
00048 
00049 #include "amesos_cholmod_internal.h"
00050 #undef ASSERT
00051 
00052 #include "metis.h"
00053 /* METIS has its own ASSERT that it reveals to the user, so remove it here: */
00054 #undef ASSERT
00055 
00056 /* and redefine it back again */
00057 #ifndef NDEBUG
00058 #define ASSERT(expression) (assert (expression))
00059 #else
00060 #define ASSERT(expression)
00061 #endif
00062 
00063 #include "amesos_cholmod_partition.h"
00064 #include "amesos_cholmod_cholesky.h"
00065 
00066 
00067 /* ========================================================================== */
00068 /* === dumpgraph ============================================================ */
00069 /* ========================================================================== */
00070 
00071 /* For dumping the input graph to METIS_NodeND, to check with METIS's onmetis
00072  * and graphchk programs.  For debugging only.  To use, uncomment this #define:
00073 #define DUMP_GRAPH
00074  */
00075 
00076 #ifdef DUMP_GRAPH
00077 #include <stdio.h>
00078 /* After dumping the graph with this routine, run "onmetis metisgraph" */
00079 static void dumpgraph (idxtype *Mp, idxtype *Mi, UF_long n,
00080     cholmod_common *Common)
00081 {
00082     UF_long i, j, p, nz ;
00083     FILE *f ;
00084     nz = Mp [n] ;
00085     printf ("Dumping METIS graph n %ld nz %ld\n", n, nz) ;    /* DUMP_GRAPH */
00086     f = fopen ("metisgraph", "w") ;
00087     if (f == NULL)
00088     {
00089   ERROR (-99, "cannot open metisgraph") ;
00090   return ;
00091     }
00092     fprintf (f, "%ld %ld\n", n, nz/2) ;         /* DUMP_GRAPH */
00093     for (j = 0 ; j < n ; j++)
00094     {
00095   for (p = Mp [j] ; p < Mp [j+1] ; p++)
00096   {
00097       i = Mi [p] ;
00098       fprintf (f, " %ld", i+1) ;          /* DUMP_GRAPH */
00099   }
00100   fprintf (f, "\n") ;           /* DUMP_GRAPH */
00101     }
00102     fclose (f) ;
00103 }
00104 #endif
00105 
00106 
00107 /* ========================================================================== */
00108 /* === metis_memory_ok ====================================================== */
00109 /* ========================================================================== */
00110 
00111 /* METIS_NodeND and METIS_NodeComputeSeparator will terminate your program it
00112  * they run out of memory.  In an attempt to workaround METIS' behavior, this
00113  * routine allocates a single block of memory of size equal to an observed
00114  * upper bound on METIS' memory usage.  It then immediately deallocates the
00115  * block.  If the allocation fails, METIS is not called.
00116  *
00117  * Median memory usage for a graph with n nodes and nz edges (counting each
00118  * edge twice, or both upper and lower triangular parts of a matrix) is
00119  * 4*nz + 40*n + 4096 integers.  A "typical" upper bound is 10*nz + 50*n + 4096
00120  * integers.  Nearly all matrices tested fit within that upper bound, with the
00121  * exception two in the UF sparse matrix collection: Schenk_IBMNA/c-64 and
00122  * Gupta/gupta2.  The latter exceeds the "upper bound" by a factor of just less
00123  * than 2.
00124  *
00125  * If you do not mind having your program terminated if it runs out of memory,
00126  * set Common->metis_memory to zero.  Its default value is 2, which allows for
00127  * some memory fragmentation, and also accounts for the Gupta/gupta2 matrix.
00128  *
00129  * An alternative, if CHOLMOD is used in MATLAB, is to use a version of METIS
00130  * (4.0.2, perhaps) proposed to George Karypis.  This version uses function
00131  * pointer for malloc and free.  They can be set to mxMalloc and mxFree
00132  * (see sputil_config in MATLAB/sputil.c).  On Linux, with gcc, you must also
00133  * compile CHOLMOD, METIS, AMD, COLAMD, and CCOLAMD with the -fexceptions
00134  * compiler flag.  With this configuration, mxMalloc safely aborts the
00135  * mexFunction, frees all memory allocted by METIS, and safely returns to
00136  * MATLAB.  You may then set Common->metis_memory = 0.
00137  */
00138 
00139 #define GUESS(nz,n) (10 * (nz) + 50 * (n) + 4096)
00140 
00141 static int metis_memory_ok
00142 (
00143     Int n,
00144     Int nz,
00145     cholmod_common *Common
00146 )
00147 {
00148     double s ;
00149     void *p ;
00150     size_t metis_guard ;
00151 
00152     if (Common->metis_memory <= 0)
00153     {
00154   /* do not prevent METIS from running out of memory */
00155   return (TRUE) ;
00156     }
00157 
00158     n  = MAX (1, n) ;
00159     nz = MAX (0, nz) ;
00160 
00161     /* compute in double, to avoid integer overflow */
00162     s = GUESS ((double) nz, (double) n) ;
00163     s *= Common->metis_memory ;
00164 
00165     if (s * sizeof (idxtype) >= ((double) Size_max))
00166     {
00167   /* don't even attempt to malloc such a large block */
00168   return (FALSE) ;
00169     }
00170 
00171     /* recompute in size_t */
00172     metis_guard = GUESS ((size_t) nz, (size_t) n) ;
00173     metis_guard *= Common->metis_memory ;
00174 
00175     /* attempt to malloc the block */
00176     p = CHOLMOD(malloc) (metis_guard, sizeof (idxtype), Common) ;
00177     if (p == NULL)
00178     {
00179   /* failure - return out-of-memory condition */
00180   return (FALSE) ;
00181     }
00182 
00183     /* success - free the block */
00184     CHOLMOD(free) (metis_guard, sizeof (idxtype), p, Common) ;
00185     return (TRUE) ;
00186 }
00187 
00188 
00189 /* ========================================================================== */
00190 /* === cholmod_metis_bisector =============================================== */
00191 /* ========================================================================== */
00192 
00193 /* Finds a set of nodes that bisects the graph of A or AA' (direct interface
00194  * to METIS_NodeComputeSeparator).
00195  *
00196  * The input matrix A must be square, symmetric (with both upper and lower
00197  * parts present) and with no diagonal entries.  These conditions are NOT
00198  * checked.
00199  */
00200 
00201 UF_long CHOLMOD(metis_bisector) /* returns separator size */
00202 (
00203     /* ---- input ---- */
00204     cholmod_sparse *A,  /* matrix to bisect */
00205     Int *Anw,   /* size A->nrow, node weights */
00206     Int *Aew,   /* size nz, edge weights */
00207     /* ---- output --- */
00208     Int *Partition, /* size A->nrow */
00209     /* --------------- */
00210     cholmod_common *Common
00211 )
00212 {
00213     Int *Ap, *Ai ;
00214     idxtype *Mp, *Mi, *Mnw, *Mew, *Mpart ;
00215     Int n, nleft, nright, j, p, csep, total_weight, lightest, nz ;
00216     int Opt [8], nn, csp ;
00217     size_t n1 ;
00218     DEBUG (Int nsep) ;
00219 
00220     /* ---------------------------------------------------------------------- */
00221     /* check inputs */
00222     /* ---------------------------------------------------------------------- */
00223 
00224     RETURN_IF_NULL_COMMON (EMPTY) ;
00225     RETURN_IF_NULL (A, EMPTY) ;
00226     RETURN_IF_NULL (Anw, EMPTY) ;
00227     RETURN_IF_NULL (Aew, EMPTY) ;
00228     RETURN_IF_NULL (Partition, EMPTY) ;
00229     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
00230     if (A->stype || A->nrow != A->ncol)
00231     {
00232   /* A must be square, with both upper and lower parts present */
00233   ERROR (CHOLMOD_INVALID, "matrix must be square, symmetric,"
00234     " and with both upper/lower parts present") ;
00235   return (EMPTY) ;
00236     }
00237     Common->status = CHOLMOD_OK ;
00238 
00239     /* ---------------------------------------------------------------------- */
00240     /* quick return */
00241     /* ---------------------------------------------------------------------- */
00242 
00243     n = A->nrow ;
00244     if (n == 0)
00245     {
00246   return (0) ;
00247     }
00248     n1 = ((size_t) n) + 1 ;
00249 
00250     /* ---------------------------------------------------------------------- */
00251     /* get inputs */
00252     /* ---------------------------------------------------------------------- */
00253 
00254     Ap = A->p ;
00255     Ai = A->i ;
00256     nz = Ap [n] ;
00257 
00258     /* ---------------------------------------------------------------------- */
00259     /* METIS does not have a 64-bit integer version */
00260     /* ---------------------------------------------------------------------- */
00261 
00262 #ifdef LONG
00263     if (sizeof (Int) > sizeof (idxtype) && MAX (n,nz) > INT_MAX / sizeof (int))
00264     {
00265   /* CHOLMOD's matrix is too large for METIS */
00266   return (EMPTY) ;
00267     }
00268 #endif
00269 
00270     /* ---------------------------------------------------------------------- */
00271     /* set default options */
00272     /* ---------------------------------------------------------------------- */
00273 
00274     Opt [0] = 0 ; /* use defaults */
00275     Opt [1] = 3 ; /* matching type */
00276     Opt [2] = 1 ; /* init. partitioning algo*/
00277     Opt [3] = 2 ; /* refinement algorithm */
00278     Opt [4] = 0 ; /* no debug */
00279     Opt [5] = 0 ; /* unused */
00280     Opt [6] = 0 ; /* unused */
00281     Opt [7] = -1 ;  /* random seed */
00282 
00283     DEBUG (for (j = 0 ; j < n ; j++) ASSERT (Anw [j] > 0)) ;
00284 
00285     /* ---------------------------------------------------------------------- */
00286     /* copy Int to METIS idxtype, if necessary */
00287     /* ---------------------------------------------------------------------- */
00288 
00289     DEBUG (for (j = 0 ; j < nz ; j++) ASSERT (Aew [j] > 0)) ;
00290     if (sizeof (Int) == sizeof (idxtype))
00291     {
00292   /* this is the typical case */
00293   Mi    = (idxtype *) Ai ;
00294   Mew   = (idxtype *) Aew ;
00295   Mp    = (idxtype *) Ap ;
00296   Mnw   = (idxtype *) Anw ;
00297   Mpart = (idxtype *) Partition ;
00298     }
00299     else
00300     {
00301   /* idxtype and Int differ; copy the graph into the METIS idxtype */
00302   Mi    = CHOLMOD(malloc) (nz, sizeof (idxtype), Common) ;
00303   Mew   = CHOLMOD(malloc) (nz, sizeof (idxtype), Common) ;
00304   Mp    = CHOLMOD(malloc) (n1, sizeof (idxtype), Common) ;
00305   Mnw   = CHOLMOD(malloc) (n,  sizeof (idxtype), Common) ;
00306   Mpart = CHOLMOD(malloc) (n,  sizeof (idxtype), Common) ;
00307   if (Common->status < CHOLMOD_OK)
00308   {
00309       CHOLMOD(free) (nz, sizeof (idxtype), Mi,    Common) ;
00310       CHOLMOD(free) (nz, sizeof (idxtype), Mew,   Common) ;
00311       CHOLMOD(free) (n1, sizeof (idxtype), Mp,    Common) ;
00312       CHOLMOD(free) (n,  sizeof (idxtype), Mnw,   Common) ;
00313       CHOLMOD(free) (n,  sizeof (idxtype), Mpart, Common) ;
00314       return (EMPTY) ;
00315   }
00316   for (p = 0 ; p < nz ; p++)
00317   {
00318       Mi [p] = Ai [p] ;
00319   }
00320   for (p = 0 ; p < nz ; p++)
00321   {
00322       Mew [p] = Aew [p] ;
00323   }
00324   for (j = 0 ; j <= n ; j++)
00325   {
00326       Mp [j] = Ap [j] ;
00327   }
00328   for (j = 0 ; j <  n ; j++)
00329   {
00330       Mnw [j] = Anw [j] ;
00331   }
00332     }
00333 
00334     /* ---------------------------------------------------------------------- */
00335     /* METIS workaround: try to ensure METIS doesn't run out of memory */
00336     /* ---------------------------------------------------------------------- */
00337 
00338     if (!metis_memory_ok (n, nz, Common))
00339     {
00340   /* METIS might ask for too much memory and thus terminate the program */
00341   if (sizeof (Int) != sizeof (idxtype))
00342   {
00343       CHOLMOD(free) (nz, sizeof (idxtype), Mi,    Common) ;
00344       CHOLMOD(free) (nz, sizeof (idxtype), Mew,   Common) ;
00345       CHOLMOD(free) (n1, sizeof (idxtype), Mp,    Common) ;
00346       CHOLMOD(free) (n,  sizeof (idxtype), Mnw,   Common) ;
00347       CHOLMOD(free) (n,  sizeof (idxtype), Mpart, Common) ;
00348   }
00349   return (EMPTY) ;
00350     }
00351 
00352     /* ---------------------------------------------------------------------- */
00353     /* partition the graph */
00354     /* ---------------------------------------------------------------------- */
00355 
00356 #ifndef NDEBUG
00357     PRINT1 (("Metis graph, n = "ID"\n", n)) ;
00358     for (j = 0 ; j < n ; j++)
00359     {
00360   Int ppp ;
00361   PRINT2 (("M(:,"ID") node weight "ID"\n", j, (Int) Mnw [j])) ;
00362   ASSERT (Mnw [j] > 0) ;
00363   for (ppp = Mp [j] ; ppp < Mp [j+1] ; ppp++)
00364   {
00365       PRINT3 ((" "ID" : "ID"\n", (Int) Mi [ppp], (Int) Mew [ppp])) ;
00366       ASSERT (Mi [ppp] != j) ;
00367       ASSERT (Mew [ppp] > 0) ;
00368   }
00369     }
00370 #endif
00371 
00372     nn = n ;
00373     METIS_NodeComputeSeparator (&nn, Mp, Mi, Mnw, Mew, Opt, &csp, Mpart) ;
00374     n = nn ;
00375     csep = csp ;
00376 
00377     PRINT1 (("METIS csep "ID"\n", csep)) ;
00378 
00379     /* ---------------------------------------------------------------------- */
00380     /* copy the results back from idxtype, if required */
00381     /* ---------------------------------------------------------------------- */
00382 
00383     if (sizeof (Int) != sizeof (idxtype))
00384     {
00385   for (j = 0 ; j < n ; j++)
00386   {
00387       Partition [j] = Mpart [j] ;
00388   }
00389   CHOLMOD(free) (nz, sizeof (idxtype), Mi,    Common) ;
00390   CHOLMOD(free) (nz, sizeof (idxtype), Mew,   Common) ;
00391   CHOLMOD(free) (n1, sizeof (idxtype), Mp,    Common) ;
00392   CHOLMOD(free) (n,  sizeof (idxtype), Mnw,   Common) ;
00393   CHOLMOD(free) (n,  sizeof (idxtype), Mpart, Common) ;
00394     }
00395 
00396     /* ---------------------------------------------------------------------- */
00397     /* ensure a reasonable separator */
00398     /* ---------------------------------------------------------------------- */
00399 
00400     /* METIS can return a valid separator with no nodes in (for example) the
00401      * left part.  In this case, there really is no separator.  CHOLMOD
00402      * prefers, in this case, for all nodes to be in the separator (and both
00403      * left and right parts to be empty).  Also, if the graph is unconnected,
00404      * METIS can return a valid empty separator.  CHOLMOD prefers at least one
00405      * node in the separator.  Note that cholmod_nested_dissection only calls
00406      * this routine on connected components, but cholmod_bisect can call this
00407      * routine for any graph. */
00408 
00409     if (csep == 0)
00410     {
00411   /* The separator is empty, select lightest node as separator.  If
00412    * ties, select the highest numbered node. */
00413   lightest = 0 ;
00414   for (j = 0 ; j < n ; j++)
00415   {
00416       if (Anw [j] <= Anw [lightest])
00417       {
00418     lightest = j ;
00419       }
00420   }
00421   PRINT1 (("Force "ID" as sep\n", lightest)) ;
00422   Partition [lightest] = 2 ;
00423   csep = Anw [lightest] ;
00424     }
00425 
00426     /* determine the node weights in the left and right part of the graph */
00427     nleft = 0 ;
00428     nright = 0 ;
00429     DEBUG (nsep = 0) ;
00430     for (j = 0 ; j < n ; j++)
00431     {
00432   PRINT1 (("Partition ["ID"] = "ID"\n", j, Partition [j])) ;
00433   if (Partition [j] == 0)
00434   {
00435       nleft += Anw [j] ;
00436   }
00437   else if (Partition [j] == 1)
00438   {
00439       nright += Anw [j] ;
00440   }
00441 #ifndef NDEBUG
00442   else
00443   {
00444       ASSERT (Partition [j] == 2) ;
00445       nsep += Anw [j] ;
00446   }
00447 #endif
00448     }
00449     ASSERT (csep == nsep) ;
00450 
00451     total_weight = nleft + nright + csep ;
00452 
00453     if (csep < total_weight)
00454     {
00455   /* The separator is less than the whole graph.  Make sure the left and
00456    * right parts are either both empty or both non-empty. */
00457   PRINT1 (("nleft "ID" nright "ID" csep "ID" tot "ID"\n",
00458     nleft, nright, csep, total_weight)) ;
00459   ASSERT (nleft + nright + csep == total_weight) ;
00460   ASSERT (nleft > 0 || nright > 0) ;
00461   if ((nleft == 0 && nright > 0) || (nleft > 0 && nright == 0))
00462   {
00463       /* left or right is empty; put all nodes in the separator */
00464       PRINT1 (("Force all in sep\n")) ;
00465       csep = total_weight ;
00466       for (j = 0 ; j < n ; j++)
00467       {
00468     Partition [j] = 2 ;
00469       }
00470   }
00471     }
00472 
00473     ASSERT (CHOLMOD(dump_partition) (n, Ap, Ai, Anw, Partition, csep, Common)) ;
00474 
00475     /* ---------------------------------------------------------------------- */
00476     /* return the sum of the weights of nodes in the separator */
00477     /* ---------------------------------------------------------------------- */
00478 
00479     return (csep) ;
00480 }
00481 
00482 
00483 /* ========================================================================== */
00484 /* === cholmod_metis ======================================================== */
00485 /* ========================================================================== */
00486 
00487 /* CHOLMOD wrapper for the METIS_NodeND ordering routine.  Creates A+A',
00488  * A*A' or A(:,f)*A(:,f)' and then calls METIS_NodeND on the resulting graph.
00489  * This routine is comparable to cholmod_nested_dissection, except that it
00490  * calls METIS_NodeND directly, and it does not return the separator tree.
00491  *
00492  * workspace:  Flag (nrow), Iwork (4*n+uncol)
00493  *  Allocates a temporary matrix B=A*A' or B=A.
00494  */
00495 
00496 int CHOLMOD(metis)
00497 (
00498     /* ---- input ---- */
00499     cholmod_sparse *A,  /* matrix to order */
00500     Int *fset,    /* subset of 0:(A->ncol)-1 */
00501     size_t fsize, /* size of fset */
00502     int postorder,  /* if TRUE, follow with etree or coletree postorder */
00503     /* ---- output --- */
00504     Int *Perm,    /* size A->nrow, output permutation */
00505     /* --------------- */
00506     cholmod_common *Common
00507 )
00508 {
00509     double d ;
00510     Int *Iperm, *Iwork, *Bp, *Bi ;
00511     idxtype *Mp, *Mi, *Mperm, *Miperm ;
00512     cholmod_sparse *B ;
00513     Int i, j, n, nz, p, identity, uncol ;
00514     int Opt [8], nn, zero = 0 ;
00515     size_t n1, s ;
00516     int ok = TRUE ;
00517 
00518     /* ---------------------------------------------------------------------- */
00519     /* get inputs */
00520     /* ---------------------------------------------------------------------- */
00521 
00522     RETURN_IF_NULL_COMMON (FALSE) ;
00523     RETURN_IF_NULL (A, FALSE) ;
00524     RETURN_IF_NULL (Perm, FALSE) ;
00525     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00526     Common->status = CHOLMOD_OK ;
00527 
00528     /* ---------------------------------------------------------------------- */
00529     /* quick return */
00530     /* ---------------------------------------------------------------------- */
00531 
00532     n = A->nrow ;
00533     if (n == 0)
00534     {
00535   return (TRUE) ;
00536     }
00537     n1 = ((size_t) n) + 1 ;
00538 
00539     /* ---------------------------------------------------------------------- */
00540     /* allocate workspace */
00541     /* ---------------------------------------------------------------------- */
00542 
00543     /* s = 4*n + uncol */
00544     uncol = (A->stype == 0) ? A->ncol : 0 ;
00545     s = CHOLMOD(mult_size_t) (n, 4, &ok) ;
00546     s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
00547     if (!ok)
00548     {
00549   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00550   return (FALSE) ;
00551     }
00552 
00553     CHOLMOD(allocate_work) (n, s, 0, Common) ;
00554     if (Common->status < CHOLMOD_OK)
00555     {
00556   return (FALSE) ;
00557     }
00558     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00559 
00560     /* ---------------------------------------------------------------------- */
00561     /* convert the matrix to adjacency list form */
00562     /* ---------------------------------------------------------------------- */
00563 
00564     /* The input graph for METIS must be symmetric, with both upper and lower
00565      * parts present, and with no diagonal entries.  The columns need not be
00566      * sorted.
00567      * B = A+A', A*A', or A(:,f)*A(:,f)', upper and lower parts present */
00568     if (A->stype)
00569     {
00570   /* Add the upper/lower part to a symmetric lower/upper matrix by
00571    * converting to unsymmetric mode */
00572   /* workspace: Iwork (nrow) */
00573   B = CHOLMOD(copy) (A, 0, -1, Common) ;
00574     }
00575     else
00576     {
00577   /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
00578   /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
00579   B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
00580     }
00581     ASSERT (CHOLMOD(dump_sparse) (B, "B for NodeND", Common) >= 0) ;
00582     if (Common->status < CHOLMOD_OK)
00583     {
00584   return (FALSE) ;
00585     }
00586     ASSERT (B->nrow == A->nrow) ;
00587 
00588     /* ---------------------------------------------------------------------- */
00589     /* get inputs */
00590     /* ---------------------------------------------------------------------- */
00591 
00592     Iwork = Common->Iwork ;
00593     Iperm = Iwork ;   /* size n (i/i/l) */
00594 
00595     Bp = B->p ;
00596     Bi = B->i ;
00597     nz = Bp [n] ;
00598 
00599     /* ---------------------------------------------------------------------- */
00600     /* METIS does not have a UF_long integer version */
00601     /* ---------------------------------------------------------------------- */
00602 
00603 #ifdef LONG
00604     if (sizeof (Int) > sizeof (idxtype) && MAX (n,nz) > INT_MAX / sizeof (int))
00605     {
00606   /* CHOLMOD's matrix is too large for METIS */
00607   CHOLMOD(free_sparse) (&B, Common) ;
00608   return (FALSE) ;
00609     }
00610 #endif
00611 
00612     /* B does not include the diagonal, and both upper and lower parts.
00613      * Common->anz includes the diagonal, and just the lower part of B */
00614     Common->anz = nz / 2 + n ;
00615 
00616     /* ---------------------------------------------------------------------- */
00617     /* set control parameters for METIS_NodeND */
00618     /* ---------------------------------------------------------------------- */
00619 
00620     Opt [0] = 0 ; /* use defaults */
00621     Opt [1] = 3 ; /* matching type */
00622     Opt [2] = 1 ; /* init. partitioning algo*/
00623     Opt [3] = 2 ; /* refinement algorithm */
00624     Opt [4] = 0 ; /* no debug */
00625     Opt [5] = 1 ; /* initial compression */
00626     Opt [6] = 0 ; /* no dense node removal */
00627     Opt [7] = 1 ; /* number of separators @ each step */
00628 
00629     /* ---------------------------------------------------------------------- */
00630     /* allocate the METIS input arrays, if needed */
00631     /* ---------------------------------------------------------------------- */
00632 
00633     if (sizeof (Int) == sizeof (idxtype))
00634     {
00635   /* This is the typical case. */
00636   Miperm = (idxtype *) Iperm ;
00637   Mperm  = (idxtype *) Perm ;
00638   Mp     = (idxtype *) Bp ;
00639   Mi     = (idxtype *) Bi ;
00640     }
00641     else
00642     {
00643   /* allocate graph for METIS only if Int and idxtype differ */
00644   Miperm = CHOLMOD(malloc) (n,  sizeof (idxtype), Common) ;
00645   Mperm  = CHOLMOD(malloc) (n,  sizeof (idxtype), Common) ;
00646   Mp     = CHOLMOD(malloc) (n1, sizeof (idxtype), Common) ;
00647   Mi     = CHOLMOD(malloc) (nz, sizeof (idxtype), Common) ;
00648   if (Common->status < CHOLMOD_OK)
00649   {
00650       /* out of memory */
00651       CHOLMOD(free_sparse) (&B, Common) ;
00652       CHOLMOD(free) (n,  sizeof (idxtype), Miperm, Common) ;
00653       CHOLMOD(free) (n,  sizeof (idxtype), Mperm, Common) ;
00654       CHOLMOD(free) (n1, sizeof (idxtype), Mp, Common) ;
00655       CHOLMOD(free) (nz, sizeof (idxtype), Mi, Common) ;
00656       return (FALSE) ;
00657   }
00658   for (j = 0 ; j <= n ; j++)
00659   {
00660       Mp [j] = Bp [j] ;
00661   }
00662   for (p = 0 ; p < nz ; p++)
00663   {
00664       Mi [p] = Bi [p] ;
00665   }
00666     }
00667 
00668     /* ---------------------------------------------------------------------- */
00669     /* METIS workarounds */
00670     /* ---------------------------------------------------------------------- */
00671 
00672     identity = FALSE ;
00673     if (nz == 0)
00674     {
00675   /* The matrix has no off-diagonal entries.  METIS_NodeND fails in this
00676    * case, so avoid using it.  The best permutation is identity anyway,
00677    * so this is an easy fix. */
00678   identity = TRUE ;
00679   PRINT1 (("METIS:: no nz\n")) ;
00680     }
00681     else if (Common->metis_nswitch > 0)
00682     {
00683   /* METIS_NodeND in METIS 4.0.1 gives a seg fault with one matrix of
00684    * order n = 3005 and nz = 6,036,025, including the diagonal entries.
00685    * The workaround is to return the identity permutation instead of using
00686    * METIS for matrices of dimension 3000 or more and with density of 66%
00687    * or more - admittedly an uncertain fix, but such matrices are so dense
00688    * that any reasonable ordering will do, even identity (n^2 is only 50%
00689    * higher than nz in this case).  CHOLMOD's nested dissection method
00690    * (cholmod_nested_dissection) has no problems with the same matrix,
00691    * even though it too uses METIS_NodeComputeSeparator.  The matrix is
00692    * derived from LPnetlib/lpi_cplex1 in the UF sparse matrix collection.
00693    * If C is the lpi_cplex matrix (of order 3005-by-5224), A = (C*C')^2
00694    * results in the seg fault.  The seg fault also occurs in the stand-
00695    * alone onmetis program that comes with METIS.  If a future version of
00696    * METIS fixes this problem, then set Common->metis_nswitch to zero.
00697    */
00698   d = ((double) nz) / (((double) n) * ((double) n)) ;
00699   if (n > (Int) (Common->metis_nswitch) && d > Common->metis_dswitch)
00700   {
00701       identity = TRUE ;
00702       PRINT1 (("METIS:: nswitch/dswitch activated\n")) ;
00703   }
00704     }
00705 
00706     if (!identity && !metis_memory_ok (n, nz, Common))
00707     {
00708   /* METIS might ask for too much memory and thus terminate the program */
00709   identity = TRUE ;
00710     }
00711 
00712     /* ---------------------------------------------------------------------- */
00713     /* find the permutation */
00714     /* ---------------------------------------------------------------------- */
00715 
00716     if (identity)
00717     {
00718   /* no need to do the postorder */
00719   postorder = FALSE ;
00720   for (i = 0 ; i < n ; i++)
00721   {
00722       Mperm [i] = i ;
00723   }
00724     }
00725     else
00726     {
00727 #ifdef DUMP_GRAPH
00728   /* DUMP_GRAPH */ printf ("Calling METIS_NodeND n "ID" nz "ID""
00729   "density %g\n", n, nz, ((double) nz) / (((double) n) * ((double) n)));
00730   dumpgraph (Mp, Mi, n, Common) ;
00731 #endif
00732 
00733   nn = n ;
00734   METIS_NodeND (&nn, Mp, Mi, &zero, Opt, Mperm, Miperm) ;
00735   n = nn ;
00736 
00737   PRINT0 (("METIS_NodeND done\n")) ;
00738     }
00739 
00740     /* ---------------------------------------------------------------------- */
00741     /* free the METIS input arrays */
00742     /* ---------------------------------------------------------------------- */
00743 
00744     if (sizeof (Int) != sizeof (idxtype))
00745     {
00746   for (i = 0 ; i < n ; i++)
00747   {
00748       Perm [i] = (Int) (Mperm [i]) ;
00749   }
00750   CHOLMOD(free) (n,   sizeof (idxtype), Miperm, Common) ;
00751   CHOLMOD(free) (n,   sizeof (idxtype), Mperm, Common) ;
00752   CHOLMOD(free) (n+1, sizeof (idxtype), Mp, Common) ;
00753   CHOLMOD(free) (nz,  sizeof (idxtype), Mi, Common) ;
00754     }
00755 
00756     CHOLMOD(free_sparse) (&B, Common) ;
00757 
00758     /* ---------------------------------------------------------------------- */
00759     /* etree or column-etree postordering, using the Cholesky Module */
00760     /* ---------------------------------------------------------------------- */
00761 
00762     if (postorder)
00763     {
00764   Int *Parent, *Post, *NewPerm ;
00765   Int k ;
00766 
00767   Parent = Iwork + 2*((size_t) n) + uncol ;   /* size n = nrow */
00768   Post   = Parent + n ;         /* size n */
00769 
00770   /* workspace: Iwork (2*nrow+uncol), Flag (nrow), Head (nrow+1) */
00771   CHOLMOD(analyze_ordering) (A, CHOLMOD_METIS, Perm, fset, fsize,
00772     Parent, Post, NULL, NULL, NULL, Common) ;
00773   if (Common->status == CHOLMOD_OK)
00774   {
00775       /* combine the METIS permutation with its postordering */
00776       NewPerm = Parent ;      /* use Parent as workspace */
00777       for (k = 0 ; k < n ; k++)
00778       {
00779     NewPerm [k] = Perm [Post [k]] ;
00780       }
00781       for (k = 0 ; k < n ; k++)
00782       {
00783     Perm [k] = NewPerm [k] ;
00784       }
00785   }
00786     }
00787 
00788     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00789     PRINT1 (("cholmod_metis done\n")) ;
00790     return (Common->status == CHOLMOD_OK) ;
00791 }
00792 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines