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