1 /* ========================================================================== */
2 /* === Partition/cholmod_metis ============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Partition Module.
7 * Copyright (C) 2005-2006, Univ. of Florida. Author: Timothy A. Davis
8 * -------------------------------------------------------------------------- */
9
10 /* CHOLMOD interface to the METIS package (Version 5.1.0):
11 *
12 * cholmod_metis_bisector:
13 *
14 * Wrapper for the METIS node separator function,
15 * METIS_ComputeVertexSeparator (METIS 5.1).
16 *
17 * Finds a set of nodes that partitions the graph into two parts. METIS
18 * 4.0 (the function METIS_ComputeVertexSeparator) allowed for edge
19 * weights to be passed to the bisector. This feature is removed in METIS
20 * 5.1. CHOLMOD itself does not rely on this feature (it calls the METIS
21 * bisector with edge weights of all 1s). However, user code can call
22 * cholmod_metis_bisector directly, and pass in edge weights. If you use
23 * METIS 5.1, these edge weights are now ignored; if you pass a non-NULL
24 * entry for edge weights, an error will be returned.
25 *
26 * cholmod_metis:
27 *
28 * Wrapper for METIS_NodeND, METIS's own nested dissection algorithm.
29 * Typically faster than cholmod_nested_dissection, mostly because it
30 * uses minimum degree on just the leaves of the separator tree, rather
31 * than the whole matrix.
32 *
33 * Note that METIS does not return an error if it runs out of memory. Instead,
34 * it terminates the program. This interface attempts to avoid that problem
35 * by preallocating space that should be large enough for any memory allocations
36 * within METIS, and then freeing that space, just before the call to METIS.
37 * While this is not guaranteed to work, it is very unlikely to fail. If you
38 * encounter this problem, increase Common->metis_memory. If you don't mind
39 * having your program terminated, set Common->metis_memory to zero (a value of
40 * 2.0 is usually safe). Several other METIS workarounds are made in the
41 * routines in this file. See the description of metis_memory_ok, just below,
42 * for more details.
43 *
44 * FUTURE WORK: interfaces to other partitioners (CHACO, SCOTCH, JOSTLE, ... )
45 *
46 * workspace: several size-nz and size-n temporary arrays. Uses no workspace
47 * in Common.
48 *
49 * Supports any xtype (pattern, real, complex, or zomplex).
50 */
51
52 #ifndef NPARTITION
53
54 #include "cholmod_internal.h"
55 #include "metis.h"
56 #include "cholmod_partition.h"
57 #include "cholmod_cholesky.h"
58
59 /* ========================================================================== */
60 /* === dumpgraph ============================================================ */
61 /* ========================================================================== */
62
63 /* For dumping the input graph to METIS_NodeND, to check with METIS's onmetis
64 * and graphchk programs. For debugging only. To use, uncomment this #define:
65 #define DUMP_GRAPH
66 */
67
68 #ifdef DUMP_GRAPH
69 #include <stdio.h>
70 /* After dumping the graph with this routine, run "onmetis metisgraph" */
dumpgraph(idx_t * Mp,idx_t * Mi,SuiteSparse_long n,cholmod_common * Common)71 static void dumpgraph (idx_t *Mp, idx_t *Mi, SuiteSparse_long n,
72 cholmod_common *Common)
73 {
74 SuiteSparse_long i, j, p, nz ;
75 FILE *f ;
76 nz = Mp [n] ;
77 printf ("Dumping METIS graph n %ld nz %ld\n", n, nz) ; /* DUMP_GRAPH */
78 f = fopen ("metisgraph", "w") ;
79 if (f == NULL)
80 {
81 ERROR (-99, "cannot open metisgraph") ;
82 return ;
83 }
84 fprintf (f, "%ld %ld\n", n, nz/2) ; /* DUMP_GRAPH */
85 for (j = 0 ; j < n ; j++)
86 {
87 for (p = Mp [j] ; p < Mp [j+1] ; p++)
88 {
89 i = Mi [p] ;
90 fprintf (f, " %ld", i+1) ; /* DUMP_GRAPH */
91 }
92 fprintf (f, "\n") ; /* DUMP_GRAPH */
93 }
94 fclose (f) ;
95 }
96 #endif
97
98
99 /* ========================================================================== */
100 /* === metis_memory_ok ====================================================== */
101 /* ========================================================================== */
102
103 /* METIS will terminate your program if
104 * they run out of memory. In an attempt to workaround METIS' behavior, this
105 * routine allocates a single block of memory of size equal to an observed
106 * upper bound on METIS' memory usage. It then immediately deallocates the
107 * block. If the allocation fails, METIS is not called.
108 *
109 * Median memory usage for a graph with n nodes and nz edges (counting each
110 * edge twice, or both upper and lower triangular parts of a matrix) is
111 * 4*nz + 40*n + 4096 integers. A "typical" upper bound is 10*nz + 50*n + 4096
112 * integers. Nearly all matrices tested fit within that upper bound, with the
113 * exception two in the UF sparse matrix collection: Schenk_IBMNA/c-64 and
114 * Gupta/gupta2. The latter exceeds the "upper bound" by a factor of just less
115 * than 2.
116 *
117 * If you do not mind having your program terminated if it runs out of memory,
118 * set Common->metis_memory to zero. Its default value is 2, which allows for
119 * some memory fragmentation, and also accounts for the Gupta/gupta2 matrix.
120 */
121
122 #define GUESS(nz,n) (10 * (nz) + 50 * (n) + 4096)
123
metis_memory_ok(Int n,Int nz,cholmod_common * Common)124 static int metis_memory_ok
125 (
126 Int n,
127 Int nz,
128 cholmod_common *Common
129 )
130 {
131 double s ;
132 void *p ;
133 size_t metis_guard ;
134
135 if (Common->metis_memory <= 0)
136 {
137 /* do not prevent METIS from running out of memory */
138 return (TRUE) ;
139 }
140
141 n = MAX (1, n) ;
142 nz = MAX (0, nz) ;
143
144 /* compute in double, to avoid integer overflow */
145 s = GUESS ((double) nz, (double) n) ;
146 s *= Common->metis_memory ;
147
148 if (s * sizeof (idx_t) >= ((double) Size_max))
149 {
150 /* don't even attempt to malloc such a large block */
151 return (FALSE) ;
152 }
153
154 /* recompute in size_t */
155 metis_guard = GUESS ((size_t) nz, (size_t) n) ;
156 metis_guard *= Common->metis_memory ;
157
158 /* attempt to malloc the block */
159 p = CHOLMOD(malloc) (metis_guard, sizeof (idx_t), Common) ;
160 if (p == NULL)
161 {
162 /* failure - return out-of-memory condition */
163 return (FALSE) ;
164 }
165
166 /* success - free the block */
167 CHOLMOD(free) (metis_guard, sizeof (idx_t), p, Common) ;
168 return (TRUE) ;
169 }
170
171
172 /* ========================================================================== */
173 /* === cholmod_metis_bisector =============================================== */
174 /* ========================================================================== */
175
176 /* Finds a set of nodes that bisects the graph of A or AA' (direct interface
177 * to METIS_ComputeVertexSeparator.
178 *
179 * The input matrix A must be square, symmetric (with both upper and lower
180 * parts present) and with no diagonal entries. These conditions are NOT
181 * checked.
182 */
183
CHOLMOD(metis_bisector)184 SuiteSparse_long CHOLMOD(metis_bisector) /* returns separator size */
185 (
186 /* ---- input ---- */
187 cholmod_sparse *A, /* matrix to bisect */
188 Int *Anw, /* size A->nrow, node weights, can be NULL, */
189 /* which means the graph is unweighted. */
190 Int *Aew, /* size nz, edge weights (silently ignored). */
191 /* This option was available with METIS 4, but not */
192 /* in METIS 5. This argument is now unused, but */
193 /* it remains for backward compatibilty, so as not */
194 /* to change the API for cholmod_metis_bisector. */
195 /* ---- output --- */
196 Int *Partition, /* size A->nrow */
197 /* --------------- */
198 cholmod_common *Common
199 )
200 {
201 Int *Ap, *Ai ;
202 idx_t *Mp, *Mi, *Mnw, *Mpart ;
203 Int n, nleft, nright, j, p, csep, total_weight, lightest, nz ;
204 idx_t nn, csp ;
205 size_t n1 ;
206 int ok ;
207 DEBUG (Int nsep) ;
208
209 /* ---------------------------------------------------------------------- */
210 /* check inputs */
211 /* ---------------------------------------------------------------------- */
212
213 RETURN_IF_NULL_COMMON (EMPTY) ;
214 RETURN_IF_NULL (A, EMPTY) ;
215 /* RETURN_IF_NULL (Anw, EMPTY) ; */
216 /* RETURN_IF_NULL (Aew, EMPTY) ; */
217 RETURN_IF_NULL (Partition, EMPTY) ;
218 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
219 if (A->stype || A->nrow != A->ncol)
220 {
221 /* A must be square, with both upper and lower parts present */
222 ERROR (CHOLMOD_INVALID, "matrix must be square, symmetric,"
223 " and with both upper/lower parts present") ;
224 return (EMPTY) ;
225 }
226 Common->status = CHOLMOD_OK ;
227 ASSERT (CHOLMOD(dump_sparse) (A, "A for bisector", Common) >= 0) ;
228
229 /* ---------------------------------------------------------------------- */
230 /* quick return */
231 /* ---------------------------------------------------------------------- */
232
233 n = A->nrow ;
234 if (n == 0)
235 {
236 return (0) ;
237 }
238 n1 = ((size_t) n) + 1 ;
239
240 /* ---------------------------------------------------------------------- */
241 /* get inputs */
242 /* ---------------------------------------------------------------------- */
243
244 Ap = A->p ;
245 Ai = A->i ;
246 nz = Ap [n] ;
247
248 if (Anw != NULL) DEBUG (for (j = 0 ; j < n ; j++) ASSERT (Anw [j] > 0)) ;
249
250 /* ---------------------------------------------------------------------- */
251 /* copy Int to METIS idx_t, if necessary */
252 /* ---------------------------------------------------------------------- */
253
254 if (sizeof (Int) == sizeof (idx_t))
255 {
256 /* this is the typical case */
257 Mi = (idx_t *) Ai ;
258 Mp = (idx_t *) Ap ;
259 Mnw = (idx_t *) Anw ;
260 Mpart = (idx_t *) Partition ;
261 }
262 else
263 {
264 /* idx_t and Int differ; copy the graph into the METIS idx_t */
265 Mi = CHOLMOD(malloc) (nz, sizeof (idx_t), Common) ;
266 Mp = CHOLMOD(malloc) (n1, sizeof (idx_t), Common) ;
267 Mnw = Anw ? (CHOLMOD(malloc) (n, sizeof (idx_t), Common)) : NULL ;
268 Mpart = CHOLMOD(malloc) (n, sizeof (idx_t), Common) ;
269 if (Common->status < CHOLMOD_OK)
270 {
271 CHOLMOD(free) (nz, sizeof (idx_t), Mi, Common) ;
272 CHOLMOD(free) (n1, sizeof (idx_t), Mp, Common) ;
273 CHOLMOD(free) (n, sizeof (idx_t), Mnw, Common) ;
274 CHOLMOD(free) (n, sizeof (idx_t), Mpart, Common) ;
275 return (EMPTY) ;
276 }
277 for (p = 0 ; p < nz ; p++)
278 {
279 Mi [p] = Ai [p] ;
280 }
281 for (j = 0 ; j <= n ; j++)
282 {
283 Mp [j] = Ap [j] ;
284 }
285 if (Anw != NULL)
286 {
287 for (j = 0 ; j < n ; j++)
288 {
289 Mnw [j] = Anw [j] ;
290 }
291 }
292 }
293
294 /* ---------------------------------------------------------------------- */
295 /* METIS workaround: try to ensure METIS doesn't run out of memory */
296 /* ---------------------------------------------------------------------- */
297
298 if (!metis_memory_ok (n, nz, Common))
299 {
300 /* METIS might ask for too much memory and thus terminate the program */
301 if (sizeof (Int) != sizeof (idx_t))
302 {
303 CHOLMOD(free) (nz, sizeof (idx_t), Mi, Common) ;
304 CHOLMOD(free) (n1, sizeof (idx_t), Mp, Common) ;
305 CHOLMOD(free) (n, sizeof (idx_t), Mnw, Common) ;
306 CHOLMOD(free) (n, sizeof (idx_t), Mpart, Common) ;
307 }
308 return (EMPTY) ;
309 }
310
311 /* ---------------------------------------------------------------------- */
312 /* partition the graph */
313 /* ---------------------------------------------------------------------- */
314
315 #ifndef NDEBUG
316 PRINT1 (("Metis graph, n = "ID"\n", n)) ;
317 for (j = 0 ; j < n ; j++)
318 {
319 Int ppp, nodeweight = (Mnw ? Mnw [j] : 1) ;
320 PRINT2 (("M(:,"ID") node weight "ID"\n", j, nodeweight)) ;
321 ASSERT (nodeweight > 0) ;
322 for (ppp = Mp [j] ; ppp < Mp [j+1] ; ppp++)
323 {
324 PRINT3 ((" "ID "\n", (Int) Mi [ppp])) ;
325 ASSERT (Mi [ppp] != j) ;
326 }
327 }
328 #endif
329
330 /*
331 METIS_ComputeVertexSeparator(
332 idx_t *nvtxs, number of nodes
333 idx_t *xadj, column pointers
334 idx_t *adjncy, row indices
335 idx_t *vwgt, vertex weights (NULL means unweighted)
336 idx_t *options, options (NULL means defaults)
337 idx_t *sepsize, separator size
338 idx_t *part); partition. part [i] = 0,1,2, where:
339 0:left, 1:right, 2:separator
340 */
341
342 nn = n ;
343 ok = METIS_ComputeVertexSeparator (&nn, Mp, Mi, Mnw, NULL, &csp, Mpart) ;
344 csep = csp ;
345
346 PRINT1 (("METIS csep "ID"\n", csep)) ;
347
348 /* ---------------------------------------------------------------------- */
349 /* copy the results back from idx_t, if required */
350 /* ---------------------------------------------------------------------- */
351
352 if (ok == METIS_OK && (sizeof (Int) != sizeof (idx_t)))
353 {
354 for (j = 0 ; j < n ; j++)
355 {
356 Partition [j] = Mpart [j] ;
357 }
358 }
359
360 /* ---------------------------------------------------------------------- */
361 /* free the workspace for METIS, if allocated */
362 /* ---------------------------------------------------------------------- */
363
364 if (sizeof (Int) != sizeof (idx_t))
365 {
366 CHOLMOD(free) (nz, sizeof (idx_t), Mi, Common) ;
367 CHOLMOD(free) (n1, sizeof (idx_t), Mp, Common) ;
368 CHOLMOD(free) (n, sizeof (idx_t), Mnw, Common) ;
369 CHOLMOD(free) (n, sizeof (idx_t), Mpart, Common) ;
370 }
371
372 if (ok == METIS_ERROR_MEMORY)
373 {
374 ERROR (CHOLMOD_OUT_OF_MEMORY, "out of memory in METIS") ;
375 return (EMPTY) ;
376 }
377 else if (ok == METIS_ERROR_INPUT)
378 {
379 ERROR (CHOLMOD_INVALID, "invalid input to METIS") ;
380 return (EMPTY) ;
381 }
382 else if (ok == METIS_ERROR)
383 {
384 ERROR (CHOLMOD_INVALID, "unspecified METIS error") ;
385 return (EMPTY) ;
386 }
387
388 /* ---------------------------------------------------------------------- */
389 /* ensure a reasonable separator */
390 /* ---------------------------------------------------------------------- */
391
392 /* METIS can return a valid separator with no nodes in (for example) the
393 * left part. In this case, there really is no separator. CHOLMOD
394 * prefers, in this case, for all nodes to be in the separator (and both
395 * left and right parts to be empty). Also, if the graph is unconnected,
396 * METIS can return a valid empty separator. CHOLMOD prefers at least one
397 * node in the separator. Note that cholmod_nested_dissection only calls
398 * this routine on connected components, but cholmod_bisect can call this
399 * routine for any graph. */
400
401 if (csep == 0)
402 {
403 /* The separator is empty, select lightest node as separator. If
404 * ties, select the highest numbered node. */
405 if (Anw == NULL)
406 {
407 lightest = n-1 ;
408 }
409 else
410 {
411 lightest = 0 ;
412 for (j = 0 ; j < n ; j++)
413 {
414 if (Anw [j] <= Anw [lightest])
415 {
416 lightest = j ;
417 }
418 }
419 }
420 PRINT1 (("Force "ID" as sep\n", lightest)) ;
421 Partition [lightest] = 2 ;
422 csep = (Anw ? (Anw [lightest]) : 1) ;
423 }
424
425 /* determine the node weights in the left and right part of the graph */
426 nleft = 0 ;
427 nright = 0 ;
428 DEBUG (nsep = 0) ;
429 for (j = 0 ; j < n ; j++)
430 {
431 PRINT1 (("Partition ["ID"] = "ID"\n", j, Partition [j])) ;
432 if (Partition [j] == 0)
433 {
434 nleft += (Anw ? (Anw [j]) : 1) ;
435 }
436 else if (Partition [j] == 1)
437 {
438 nright += (Anw ? (Anw [j]) : 1) ;
439 }
440 #ifndef NDEBUG
441 else
442 {
443 ASSERT (Partition [j] == 2) ;
444 nsep += (Anw ? (Anw [j]) : 1) ;
445 }
446 #endif
447 }
448 ASSERT (csep == nsep) ;
449
450 total_weight = nleft + nright + csep ;
451
452 if (csep < total_weight)
453 {
454 /* The separator is less than the whole graph. Make sure the left and
455 * right parts are either both empty or both non-empty. */
456 PRINT1 (("nleft "ID" nright "ID" csep "ID" tot "ID"\n",
457 nleft, nright, csep, total_weight)) ;
458 ASSERT (nleft + nright + csep == total_weight) ;
459 ASSERT (nleft > 0 || nright > 0) ;
460 if ((nleft == 0 && nright > 0) || (nleft > 0 && nright == 0))
461 {
462 /* left or right is empty; put all nodes in the separator */
463 PRINT1 (("Force all in sep\n")) ;
464 csep = total_weight ;
465 for (j = 0 ; j < n ; j++)
466 {
467 Partition [j] = 2 ;
468 }
469 }
470 }
471
472 ASSERT (CHOLMOD(dump_partition) (n, Ap, Ai, Anw, Partition, csep, Common)) ;
473
474 /* ---------------------------------------------------------------------- */
475 /* return the sum of the weights of nodes in the separator */
476 /* ---------------------------------------------------------------------- */
477
478 return (csep) ;
479 }
480
481
482 /* ========================================================================== */
483 /* === cholmod_metis ======================================================== */
484 /* ========================================================================== */
485
486 /* CHOLMOD wrapper for the METIS_NodeND ordering routine. Creates A+A',
487 * A*A' or A(:,f)*A(:,f)' and then calls METIS_NodeND on the resulting graph.
488 * This routine is comparable to cholmod_nested_dissection, except that it
489 * calls METIS_NodeND directly, and it does not return the separator tree.
490 *
491 * workspace: Flag (nrow), Iwork (4*n+uncol)
492 * Allocates a temporary matrix B=A*A' or B=A.
493 */
494
CHOLMOD(metis)495 int CHOLMOD(metis)
496 (
497 /* ---- input ---- */
498 cholmod_sparse *A, /* matrix to order */
499 Int *fset, /* subset of 0:(A->ncol)-1 */
500 size_t fsize, /* size of fset */
501 int postorder, /* if TRUE, follow with etree or coletree postorder */
502 /* ---- output --- */
503 Int *Perm, /* size A->nrow, output permutation */
504 /* --------------- */
505 cholmod_common *Common
506 )
507 {
508 double d ;
509 Int *Iperm, *Iwork, *Bp, *Bi ;
510 idx_t *Mp, *Mi, *Mperm, *Miperm ;
511 cholmod_sparse *B ;
512 Int i, j, n, nz, p, identity, uncol ;
513 idx_t nn, zero = 0 ;
514 size_t n1, s ;
515 int ok = TRUE ;
516
517 /* ---------------------------------------------------------------------- */
518 /* get inputs */
519 /* ---------------------------------------------------------------------- */
520
521 RETURN_IF_NULL_COMMON (FALSE) ;
522 RETURN_IF_NULL (A, FALSE) ;
523 RETURN_IF_NULL (Perm, FALSE) ;
524 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
525 Common->status = CHOLMOD_OK ;
526
527 /* ---------------------------------------------------------------------- */
528 /* quick return */
529 /* ---------------------------------------------------------------------- */
530
531 n = A->nrow ;
532 if (n == 0)
533 {
534 return (TRUE) ;
535 }
536 n1 = ((size_t) n) + 1 ;
537
538 /* ---------------------------------------------------------------------- */
539 /* allocate workspace */
540 /* ---------------------------------------------------------------------- */
541
542 /* s = 4*n + uncol */
543 uncol = (A->stype == 0) ? A->ncol : 0 ;
544 s = CHOLMOD(mult_size_t) (n, 4, &ok) ;
545 s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
546 if (!ok)
547 {
548 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
549 return (FALSE) ;
550 }
551
552 CHOLMOD(allocate_work) (n, s, 0, Common) ;
553 if (Common->status < CHOLMOD_OK)
554 {
555 return (FALSE) ;
556 }
557 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
558
559 /* ---------------------------------------------------------------------- */
560 /* convert the matrix to adjacency list form */
561 /* ---------------------------------------------------------------------- */
562
563 /* The input graph for METIS must be symmetric, with both upper and lower
564 * parts present, and with no diagonal entries. The columns need not be
565 * sorted.
566 * B = A+A', A*A', or A(:,f)*A(:,f)', upper and lower parts present */
567 if (A->stype)
568 {
569 /* Add the upper/lower part to a symmetric lower/upper matrix by
570 * converting to unsymmetric mode */
571 /* workspace: Iwork (nrow) */
572 B = CHOLMOD(copy) (A, 0, -1, Common) ;
573 }
574 else
575 {
576 /* B = A*A' or A(:,f)*A(:,f)', no diagonal */
577 /* workspace: Flag (nrow), Iwork (max (nrow,ncol)) */
578 B = CHOLMOD(aat) (A, fset, fsize, -1, Common) ;
579 }
580 ASSERT (CHOLMOD(dump_sparse) (B, "B for NodeND", Common) >= 0) ;
581 if (Common->status < CHOLMOD_OK)
582 {
583 return (FALSE) ;
584 }
585 ASSERT (B->nrow == A->nrow) ;
586
587 /* ---------------------------------------------------------------------- */
588 /* get inputs */
589 /* ---------------------------------------------------------------------- */
590
591 Iwork = Common->Iwork ;
592 Iperm = Iwork ; /* size n (i/i/l) */
593
594 Bp = B->p ;
595 Bi = B->i ;
596 nz = Bp [n] ;
597
598 /* B does not include the diagonal, and both upper and lower parts.
599 * Common->anz includes the diagonal, and just the lower part of B */
600 Common->anz = nz / 2 + n ;
601
602 /* ---------------------------------------------------------------------- */
603 /* allocate the METIS input arrays, if needed */
604 /* ---------------------------------------------------------------------- */
605
606 if (sizeof (Int) == sizeof (idx_t))
607 {
608 /* This is the typical case. */
609 Miperm = (idx_t *) Iperm ;
610 Mperm = (idx_t *) Perm ;
611 Mp = (idx_t *) Bp ;
612 Mi = (idx_t *) Bi ;
613 }
614 else
615 {
616 /* allocate graph for METIS only if Int and idx_t differ */
617 Miperm = CHOLMOD(malloc) (n, sizeof (idx_t), Common) ;
618 Mperm = CHOLMOD(malloc) (n, sizeof (idx_t), Common) ;
619 Mp = CHOLMOD(malloc) (n1, sizeof (idx_t), Common) ;
620 Mi = CHOLMOD(malloc) (nz, sizeof (idx_t), Common) ;
621 if (Common->status < CHOLMOD_OK)
622 {
623 /* out of memory */
624 CHOLMOD(free_sparse) (&B, Common) ;
625 CHOLMOD(free) (n, sizeof (idx_t), Miperm, Common) ;
626 CHOLMOD(free) (n, sizeof (idx_t), Mperm, Common) ;
627 CHOLMOD(free) (n1, sizeof (idx_t), Mp, Common) ;
628 CHOLMOD(free) (nz, sizeof (idx_t), Mi, Common) ;
629 return (FALSE) ;
630 }
631 for (j = 0 ; j <= n ; j++)
632 {
633 Mp [j] = Bp [j] ;
634 }
635 for (p = 0 ; p < nz ; p++)
636 {
637 Mi [p] = Bi [p] ;
638 }
639 }
640
641 /* ---------------------------------------------------------------------- */
642 /* METIS workarounds */
643 /* ---------------------------------------------------------------------- */
644
645 identity = FALSE ;
646 if (nz == 0)
647 {
648 /* The matrix has no off-diagonal entries. METIS_NodeND fails in this
649 * case, so avoid using it. The best permutation is identity anyway,
650 * so this is an easy fix. */
651 identity = TRUE ;
652 PRINT1 (("METIS:: no nz\n")) ;
653 }
654 else if (Common->metis_nswitch > 0)
655 {
656 /* METIS_NodeND in METIS 4.0.1 gives a seg fault with one matrix of
657 * order n = 3005 and nz = 6,036,025, including the diagonal entries.
658 * The workaround is to return the identity permutation instead of using
659 * METIS for matrices of dimension 3000 or more and with density of 66%
660 * or more - admittedly an uncertain fix, but such matrices are so dense
661 * that any reasonable ordering will do, even identity (n^2 is only 50%
662 * higher than nz in this case). CHOLMOD's nested dissection method
663 * (cholmod_nested_dissection) has no problems with the same matrix,
664 * even though it too uses METIS_ComputeVertexSeparator. The matrix is
665 * derived from LPnetlib/lpi_cplex1 in the UF sparse matrix collection.
666 * If C is the lpi_cplex matrix (of order 3005-by-5224), A = (C*C')^2
667 * results in the seg fault. The seg fault also occurs in the stand-
668 * alone onmetis program that comes with METIS. If a future version of
669 * METIS fixes this problem, then set Common->metis_nswitch to zero.
670 */
671 d = ((double) nz) / (((double) n) * ((double) n)) ;
672 if (n > (Int) (Common->metis_nswitch) && d > Common->metis_dswitch)
673 {
674 identity = TRUE ;
675 PRINT1 (("METIS:: nswitch/dswitch activated\n")) ;
676 }
677 }
678
679 if (!identity && !metis_memory_ok (n, nz, Common))
680 {
681 /* METIS might ask for too much memory and thus terminate the program */
682 identity = TRUE ;
683 }
684
685 /* ---------------------------------------------------------------------- */
686 /* find the permutation */
687 /* ---------------------------------------------------------------------- */
688
689 if (identity)
690 {
691 /* no need to do the postorder */
692 postorder = FALSE ;
693 for (i = 0 ; i < n ; i++)
694 {
695 Mperm [i] = i ;
696 }
697 }
698 else
699 {
700 #ifdef DUMP_GRAPH
701 /* DUMP_GRAPH */ printf ("Calling METIS_NodeND n "ID" nz "ID""
702 "density %g\n", n, nz, ((double) nz) / (((double) n) * ((double) n)));
703 dumpgraph (Mp, Mi, n, Common) ;
704 #endif
705
706 /*
707 int METIS_NodeND(
708 idx_t *nvtxs, number of nodes
709 idx_t *xadj, column pointers
710 idx_t *adjncy, row indices
711 idx_t *vwgt, vertex weights (NULL means unweighted)
712 idx_t *options, options (NULL means defaults)
713 idx_t *perm, fill-reducing ordering
714 idx_t *iperm); inverse of perm
715 */
716
717 nn = n ;
718 METIS_NodeND (&nn, Mp, Mi, NULL, NULL, Mperm, Miperm) ;
719
720 PRINT0 (("METIS_NodeND done\n")) ;
721 }
722
723 /* ---------------------------------------------------------------------- */
724 /* free the METIS input arrays */
725 /* ---------------------------------------------------------------------- */
726
727 if (sizeof (Int) != sizeof (idx_t))
728 {
729 for (i = 0 ; i < n ; i++)
730 {
731 Perm [i] = (Int) (Mperm [i]) ;
732 }
733 CHOLMOD(free) (n, sizeof (idx_t), Miperm, Common) ;
734 CHOLMOD(free) (n, sizeof (idx_t), Mperm, Common) ;
735 CHOLMOD(free) (n+1, sizeof (idx_t), Mp, Common) ;
736 CHOLMOD(free) (nz, sizeof (idx_t), Mi, Common) ;
737 }
738
739 CHOLMOD(free_sparse) (&B, Common) ;
740
741 /* ---------------------------------------------------------------------- */
742 /* etree or column-etree postordering, using the Cholesky Module */
743 /* ---------------------------------------------------------------------- */
744
745 if (postorder)
746 {
747 Int *Parent, *Post, *NewPerm ;
748 Int k ;
749
750 Parent = Iwork + 2*((size_t) n) + uncol ; /* size n = nrow */
751 Post = Parent + n ; /* size n */
752
753 /* workspace: Iwork (2*nrow+uncol), Flag (nrow), Head (nrow+1) */
754 CHOLMOD(analyze_ordering) (A, CHOLMOD_METIS, Perm, fset, fsize,
755 Parent, Post, NULL, NULL, NULL, Common) ;
756 if (Common->status == CHOLMOD_OK)
757 {
758 /* combine the METIS permutation with its postordering */
759 NewPerm = Parent ; /* use Parent as workspace */
760 for (k = 0 ; k < n ; k++)
761 {
762 NewPerm [k] = Perm [Post [k]] ;
763 }
764 for (k = 0 ; k < n ; k++)
765 {
766 Perm [k] = NewPerm [k] ;
767 }
768 }
769 }
770
771 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
772 PRINT1 (("cholmod_metis done\n")) ;
773 return (Common->status == CHOLMOD_OK) ;
774 }
775 #endif
776