1 /* ========================================================================== */
2 /* === Cholesky/cholmod_analyze ============================================= */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2013, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Order and analyze a matrix (either simplicial or supernodal), in prepartion
10 * for numerical factorization via cholmod_factorize or via the "expert"
11 * routines cholmod_rowfac and cholmod_super_numeric.
12 *
13 * symmetric case: A or A(p,p)
14 * unsymmetric case: AA', A(p,:)*A(p,:)', A(:,f)*A(:,f)', or A(p,f)*A(p,f)'
15 *
16 * For the symmetric case, only the upper or lower triangular part of A is
17 * accessed (depending on the type of A). LL'=A (or permuted A) is analzed.
18 * For the unsymmetric case (LL'=AA' or permuted A).
19 *
20 * There can be no duplicate entries in p or f. p is of length m if A is
21 * m-by-n. f can be length 0 to n.
22 *
23 * In both cases, the columns of A need not be sorted. A can be in packed
24 * or unpacked form.
25 *
26 * Ordering options include:
27 *
28 * natural: A is not permuted to reduce fill-in
29 * given: a permutation can be provided to this routine (UserPerm)
30 * AMD: approximate minumum degree (AMD for the symmetric case,
31 * COLAMD for the AA' case).
32 * METIS: nested dissection with METIS_NodeND
33 * NESDIS: nested dissection using METIS_ComputeVertexSeparator,
34 * typically followed by a constrained minimum degree
35 * (CAMD for the symmetric case, CCOLAMD for the AA' case).
36 *
37 * Multiple ordering options can be tried (up to 9 of them), and the best one
38 * is selected (the one that gives the smallest number of nonzeros in the
39 * simplicial factor L). If one method fails, cholmod_analyze keeps going, and
40 * picks the best among the methods that succeeded. This routine fails (and
41 * returns NULL) if either initial memory allocation fails, all ordering methods
42 * fail, or the supernodal analysis (if requested) fails. By default, the 9
43 * methods available are:
44 *
45 * 1) given permutation (skipped if UserPerm is NULL)
46 * 2) AMD (symmetric case) or COLAMD (unsymmetric case)
47 * 3) METIS with default parameters
48 * 4) NESDIS with default parameters (stopping the partitioning when
49 * the graph is of size nd_small = 200 or less, remove nodes with
50 * more than max (16, prune_dense * sqrt (n)) nodes where
51 * prune_dense = 10, and follow partitioning with CCOLAMD, a
52 * constrained minimum degree ordering).
53 * 5) natural
54 * 6) NESDIS, nd_small = 20000, prune_dense = 10
55 * 7) NESDIS, nd_small = 4, prune_dense = 10, no min degree
56 * 8) NESDIS, nd_small = 200, prune_dense = 0
57 * 9) COLAMD for A*A' or AMD for A
58 *
59 * By default, the first two are tried, and METIS is tried if AMD reports a high
60 * flop count and fill-in. Let fl denote the flop count for the AMD, ordering,
61 * nnz(L) the # of nonzeros in L, and nnz(tril(A)) (or A*A'). If
62 * fl/nnz(L) >= 500 and nnz(L)/nnz(tril(A)) >= 5, then METIS is attempted. The
63 * best ordering is used (UserPerm if given, AMD, and METIS if attempted). If
64 * you do not have METIS, only the first two will be tried (user permutation,
65 * if provided, and AMD/COLAMD). This default behavior is obtained when
66 * Common->nmethods is zero. In this case, methods 0, 1, and 2 in
67 * Common->method [..] are reset to User-provided, AMD, and METIS (or NESDIS
68 * if Common->default_nesdis is set to the non-default value of TRUE),
69 * respectively.
70 *
71 * You can modify these 9 methods and the number of methods tried by changing
72 * parameters in the Common argument. If you know the best ordering for your
73 * matrix, set Common->nmethods to 1 and set Common->method[0].ordering to the
74 * requested ordering method. Parameters for each method can also be modified
75 * (refer to cholmod.h for details).
76 *
77 * Note that it is possible for METIS to terminate your program if it runs out
78 * of memory. This is not the case for any CHOLMOD or minimum degree ordering
79 * routine (AMD, COLAMD, CAMD, CCOLAMD, or CSYMAMD). Since NESDIS relies on
80 * METIS, it too can terminate your program.
81 *
82 * The factor L is returned as simplicial symbolic (L->is_super FALSE) if
83 * Common->supernodal <= CHOLMOD_SIMPLICIAL (0) or as supernodal symbolic if
84 * Common->supernodal >= CHOLMOD_SUPERNODAL (2). If Common->supernodal is
85 * equal to CHOLMOD_AUTO (1), then L is simplicial if the flop count per
86 * nonzero in L is less than Common->supernodal_switch (default: 40), and
87 * is returned as a supernodal factor otherwise.
88 *
89 * In both cases, L->xtype is CHOLMOD_PATTERN.
90 * A subsequent call to cholmod_factorize will perform a
91 * simplicial or supernodal factorization, depending on the type of L.
92 *
93 * For the simplicial case, L contains the fill-reducing permutation (L->Perm)
94 * and the counts of nonzeros in each column of L (L->ColCount). For the
95 * supernodal case, L also contains the nonzero pattern of each supernode.
96 *
97 * workspace: Flag (nrow), Head (nrow+1)
98 * if symmetric: Iwork (6*nrow)
99 * if unsymmetric: Iwork (6*nrow+ncol).
100 * calls various ordering routines, which typically allocate O(nnz(A))
101 * temporary workspace ((2 to 3)*nnz(A) * sizeof (Int) is typical, but it
102 * can be much higher if A*A' must be explicitly formed for METIS). Also
103 * allocates up to 2 temporary (permuted/transpose) copies of the nonzero
104 * pattern of A, and up to 3*n*sizeof(Int) additional workspace.
105 *
106 * Supports any xtype (pattern, real, complex, or zomplex)
107 */
108
109 #ifndef NCHOLESKY
110
111 #include "cholmod_internal.h"
112 #include "cholmod_cholesky.h"
113
114 #ifndef NSUPERNODAL
115 #include "cholmod_supernodal.h"
116 #endif
117
118 #ifndef NPARTITION
119 #include "cholmod_partition.h"
120 #endif
121
122
123 /* ========================================================================== */
124 /* === cholmod_analyze ====================================================== */
125 /* ========================================================================== */
126
127 /* Orders and analyzes A, AA', PAP', or PAA'P' and returns a symbolic factor
128 * that can later be passed to cholmod_factorize. */
129
CHOLMOD(analyze)130 cholmod_factor *CHOLMOD(analyze)
131 (
132 /* ---- input ---- */
133 cholmod_sparse *A, /* matrix to order and analyze */
134 /* --------------- */
135 cholmod_common *Common
136 )
137 {
138 return (CHOLMOD(analyze_p2) (TRUE, A, NULL, NULL, 0, Common)) ;
139 }
140
141
142 /* ========================================================================== */
143 /* === cholmod_analyze_p ==================================================== */
144 /* ========================================================================== */
145
146 /* Orders and analyzes A, AA', PAP', PAA'P', FF', or PFF'P and returns a
147 * symbolic factor that can later be passed to cholmod_factorize, where
148 * F = A(:,fset) if fset is not NULL and A->stype is zero.
149 * UserPerm is tried if non-NULL. */
150
CHOLMOD(analyze_p)151 cholmod_factor *CHOLMOD(analyze_p)
152 (
153 /* ---- input ---- */
154 cholmod_sparse *A, /* matrix to order and analyze */
155 Int *UserPerm, /* user-provided permutation, size A->nrow */
156 Int *fset, /* subset of 0:(A->ncol)-1 */
157 size_t fsize, /* size of fset */
158 /* --------------- */
159 cholmod_common *Common
160 )
161 {
162 return (CHOLMOD(analyze_p2) (TRUE, A, UserPerm, fset, fsize, Common)) ;
163 }
164
165
166 /* ========================================================================== */
167 /* === permute_matrices ===================================================== */
168 /* ========================================================================== */
169
170 /* Permute and transpose a matrix. Allocates the A1 and A2 matrices, if needed,
171 * or returns them as NULL if not needed.
172 */
173
permute_matrices(cholmod_sparse * A,Int ordering,Int * Perm,Int * fset,size_t fsize,Int do_rowcolcounts,cholmod_sparse ** A1_handle,cholmod_sparse ** A2_handle,cholmod_sparse ** S_handle,cholmod_sparse ** F_handle,cholmod_common * Common)174 static int permute_matrices
175 (
176 /* ---- input ---- */
177 cholmod_sparse *A, /* matrix to permute */
178 Int ordering, /* ordering method used */
179 Int *Perm, /* fill-reducing permutation */
180 Int *fset, /* subset of 0:(A->ncol)-1 */
181 size_t fsize, /* size of fset */
182 Int do_rowcolcounts,/* if TRUE, compute both S and F. If FALSE, only
183 * S is needed for the symmetric case, and only F for
184 * the unsymmetric case */
185 /* ---- output --- */
186 cholmod_sparse **A1_handle, /* see comments below for A1, A2, S, F */
187 cholmod_sparse **A2_handle,
188 cholmod_sparse **S_handle,
189 cholmod_sparse **F_handle,
190 /* --------------- */
191 cholmod_common *Common
192 )
193 {
194 cholmod_sparse *A1, *A2, *S, *F ;
195
196 *A1_handle = NULL ;
197 *A2_handle = NULL ;
198 *S_handle = NULL ;
199 *F_handle = NULL ;
200 A1 = NULL ;
201 A2 = NULL ;
202
203 if (ordering == CHOLMOD_NATURAL)
204 {
205
206 /* ------------------------------------------------------------------ */
207 /* natural ordering of A */
208 /* ------------------------------------------------------------------ */
209
210 if (A->stype < 0)
211 {
212 /* symmetric lower case: A already in lower form, so S=A' */
213 /* workspace: Iwork (nrow) */
214 A2 = CHOLMOD(ptranspose) (A, 0, NULL, NULL, 0, Common) ;
215 F = A ;
216 S = A2 ;
217 }
218 else if (A->stype > 0)
219 {
220 /* symmetric upper case: F = pattern of triu (A)', S = A */
221 /* workspace: Iwork (nrow) */
222 if (do_rowcolcounts)
223 {
224 /* F not needed for symmetric case if do_rowcolcounts FALSE */
225 A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
226 }
227 F = A1 ;
228 S = A ;
229 }
230 else
231 {
232 /* unsymmetric case: F = pattern of A (:,f)', S = A */
233 /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
234 A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
235 F = A1 ;
236 S = A ;
237 }
238
239 }
240 else
241 {
242
243 /* ------------------------------------------------------------------ */
244 /* A is permuted */
245 /* ------------------------------------------------------------------ */
246
247 if (A->stype < 0)
248 {
249 /* symmetric lower case: S = tril (A (p,p))' and F = S' */
250 /* workspace: Iwork (2*nrow) */
251 A2 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
252 S = A2 ;
253 /* workspace: Iwork (nrow) */
254 if (do_rowcolcounts)
255 {
256 /* F not needed for symmetric case if do_rowcolcounts FALSE */
257 A1 = CHOLMOD(ptranspose) (A2, 0, NULL, NULL, 0, Common) ;
258 }
259 F = A1 ;
260 }
261 else if (A->stype > 0)
262 {
263 /* symmetric upper case: F = triu (A (p,p))' and S = F' */
264 /* workspace: Iwork (2*nrow) */
265 A1 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
266 F = A1 ;
267 /* workspace: Iwork (nrow) */
268 A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
269 S = A2 ;
270 }
271 else
272 {
273 /* unsymmetric case: F = A (p,f)' and S = F' */
274 /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
275 A1 = CHOLMOD(ptranspose) (A, 0, Perm, fset, fsize, Common) ;
276 F = A1 ;
277 if (do_rowcolcounts)
278 {
279 /* S not needed for unsymmetric case if do_rowcolcounts FALSE */
280 /* workspace: Iwork (nrow) */
281 A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
282 }
283 S = A2 ;
284 }
285 }
286
287 /* If any cholmod_*transpose fails, one or more matrices will be NULL */
288 *A1_handle = A1 ;
289 *A2_handle = A2 ;
290 *S_handle = S ;
291 *F_handle = F ;
292 return (Common->status == CHOLMOD_OK) ;
293 }
294
295
296 /* ========================================================================== */
297 /* === cholmod_analyze_ordering ============================================= */
298 /* ========================================================================== */
299
300 /* Given a matrix A and its fill-reducing permutation, compute the elimination
301 * tree, its (non-weighted) postordering, and the number of nonzeros in each
302 * column of L. Also computes the flop count, the total nonzeros in L, and
303 * the nonzeros in A (Common->fl, Common->lnz, and Common->anz).
304 *
305 * The column counts of L, flop count, and other statistics from
306 * cholmod_rowcolcounts are not computed if ColCount is NULL.
307 *
308 * workspace: Iwork (2*nrow if symmetric, 2*nrow+ncol if unsymmetric),
309 * Flag (nrow), Head (nrow+1)
310 */
311
CHOLMOD(analyze_ordering)312 int CHOLMOD(analyze_ordering)
313 (
314 /* ---- input ---- */
315 cholmod_sparse *A, /* matrix to analyze */
316 int ordering, /* ordering method used */
317 Int *Perm, /* size n, fill-reducing permutation to analyze */
318 Int *fset, /* subset of 0:(A->ncol)-1 */
319 size_t fsize, /* size of fset */
320 /* ---- output --- */
321 Int *Parent, /* size n, elimination tree */
322 Int *Post, /* size n, postordering of elimination tree */
323 Int *ColCount, /* size n, nnz in each column of L */
324 /* ---- workspace */
325 Int *First, /* size n workspace for cholmod_postorder */
326 Int *Level, /* size n workspace for cholmod_postorder */
327 /* --------------- */
328 cholmod_common *Common
329 )
330 {
331 cholmod_sparse *A1, *A2, *S, *F ;
332 Int n, ok, do_rowcolcounts ;
333
334 /* check inputs */
335 RETURN_IF_NULL_COMMON (FALSE) ;
336 RETURN_IF_NULL (A, FALSE) ;
337
338 n = A->nrow ;
339
340 do_rowcolcounts = (ColCount != NULL) ;
341
342 /* permute A according to Perm and fset */
343 ok = permute_matrices (A, ordering, Perm, fset, fsize, do_rowcolcounts,
344 &A1, &A2, &S, &F, Common) ;
345
346 /* find etree of S (symmetric upper/lower case) or F (unsym case) */
347 /* workspace: symmmetric: Iwork (nrow), unsym: Iwork (nrow+ncol) */
348 ok = ok && CHOLMOD(etree) (A->stype ? S:F, Parent, Common) ;
349
350 /* postorder the etree (required by cholmod_rowcolcounts) */
351 /* workspace: Iwork (2*nrow) */
352 ok = ok && (CHOLMOD(postorder) (Parent, n, NULL, Post, Common) == n) ;
353
354 /* cholmod_postorder doesn't set Common->status if it returns < n */
355 Common->status = (!ok && Common->status == CHOLMOD_OK) ?
356 CHOLMOD_INVALID : Common->status ;
357
358 /* analyze LL'=S or SS' or S(:,f)*S(:,f)' */
359 /* workspace:
360 * if symmetric: Flag (nrow), Iwork (2*nrow)
361 * if unsymmetric: Flag (nrow), Iwork (2*nrow+ncol), Head (nrow+1)
362 */
363 if (do_rowcolcounts)
364 {
365 ok = ok && CHOLMOD(rowcolcounts) (A->stype ? F:S, fset, fsize, Parent,
366 Post, NULL, ColCount, First, Level, Common) ;
367 }
368
369 /* free temporary matrices and return result */
370 CHOLMOD(free_sparse) (&A1, Common) ;
371 CHOLMOD(free_sparse) (&A2, Common) ;
372 return (ok) ;
373 }
374
375
376 /* ========================================================================== */
377 /* === Free workspace and return L ========================================== */
378 /* ========================================================================== */
379
380 #define FREE_WORKSPACE_AND_RETURN \
381 { \
382 Common->no_workspace_reallocate = FALSE ; \
383 CHOLMOD(free) (n, sizeof (Int), Lparent, Common) ; \
384 CHOLMOD(free) (n, sizeof (Int), Perm, Common) ; \
385 CHOLMOD(free) (n, sizeof (Int), ColCount, Common) ; \
386 if (Common->status < CHOLMOD_OK) \
387 { \
388 CHOLMOD(free_factor) (&L, Common) ; \
389 } \
390 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ; \
391 return (L) ; \
392 }
393
394
395 /* ========================================================================== */
396 /* === cholmod_analyze_p2 =================================================== */
397 /* ========================================================================== */
398
399 /* Ordering and analysis for sparse Cholesky or sparse QR. */
400
CHOLMOD(analyze_p2)401 cholmod_factor *CHOLMOD(analyze_p2)
402 (
403 /* ---- input ---- */
404 int for_whom, /* FOR_SPQR (0): for SPQR but not GPU-accelerated
405 FOR_CHOLESKY (1): for Cholesky (GPU or not)
406 FOR_SPQRGPU (2): for SPQR with GPU acceleration */
407 cholmod_sparse *A, /* matrix to order and analyze */
408 Int *UserPerm, /* user-provided permutation, size A->nrow */
409 Int *fset, /* subset of 0:(A->ncol)-1 */
410 size_t fsize, /* size of fset */
411 /* --------------- */
412 cholmod_common *Common
413 )
414 {
415 double lnz_best ;
416 Int *First, *Level, *Work4n, *Cmember, *CParent, *ColCount, *Lperm, *Parent,
417 *Post, *Perm, *Lparent, *Lcolcount ;
418 cholmod_factor *L ;
419 Int k, n, ordering, method, nmethods, status, default_strategy, ncol, uncol,
420 skip_analysis, skip_best ;
421 Int amd_backup ;
422 size_t s ;
423 int ok = TRUE ;
424
425 /* ---------------------------------------------------------------------- */
426 /* check inputs */
427 /* ---------------------------------------------------------------------- */
428
429 RETURN_IF_NULL_COMMON (NULL) ;
430 RETURN_IF_NULL (A, NULL) ;
431 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
432 Common->status = CHOLMOD_OK ;
433 status = CHOLMOD_OK ;
434 Common->selected = EMPTY ;
435 Common->called_nd = FALSE ;
436
437 /* ---------------------------------------------------------------------- */
438 /* get inputs */
439 /* ---------------------------------------------------------------------- */
440
441 n = A->nrow ;
442 ncol = A->ncol ;
443 uncol = (A->stype == 0) ? (A->ncol) : 0 ;
444
445 /* ---------------------------------------------------------------------- */
446 /* set the default strategy */
447 /* ---------------------------------------------------------------------- */
448
449 lnz_best = (double) EMPTY ;
450 skip_best = FALSE ;
451 nmethods = MIN (Common->nmethods, CHOLMOD_MAXMETHODS) ;
452 nmethods = MAX (0, nmethods) ;
453
454 #ifndef NDEBUG
455 PRINT1 (("cholmod_analyze_p2 :: nmethods "ID"\n", nmethods)) ;
456 for (method = 0 ; method < nmethods ; method++)
457 {
458 PRINT1 ((" "ID": ordering "ID"\n",
459 method, Common->method [method].ordering)) ;
460 }
461 #endif
462
463 default_strategy = (nmethods == 0) ;
464 if (default_strategy)
465 {
466 /* default strategy: try UserPerm, if given. Try AMD for A, or AMD
467 * to order A*A'. Try METIS for the symmetric case only if AMD reports
468 * a high degree of fill-in and flop count. METIS is not tried if the
469 * Partition Module isn't installed. If Common->default_nesdis is
470 * TRUE, then NESDIS is used as the 3rd ordering instead. */
471 Common->method [0].ordering = CHOLMOD_GIVEN ;/* skip if UserPerm NULL */
472 Common->method [1].ordering = CHOLMOD_AMD ;
473 Common->method [2].ordering =
474 (Common->default_nesdis ? CHOLMOD_NESDIS : CHOLMOD_METIS) ;
475 amd_backup = FALSE ;
476 #ifndef NPARTITION
477 nmethods = 3 ;
478 #else
479 nmethods = 2 ;
480 #endif
481 }
482 else
483 {
484 /* If only METIS and NESDIS are selected, or if 2 or more methods are
485 * being tried, then enable AMD backup */
486 amd_backup = (nmethods > 1) || (nmethods == 1 &&
487 (Common->method [0].ordering == CHOLMOD_METIS ||
488 Common->method [0].ordering == CHOLMOD_NESDIS)) ;
489 }
490
491 #ifdef NSUPERNODAL
492 /* CHOLMOD Supernodal module not installed, just do simplicial analysis */
493 Common->supernodal = CHOLMOD_SIMPLICIAL ;
494 #endif
495
496 /* ---------------------------------------------------------------------- */
497 /* allocate workspace */
498 /* ---------------------------------------------------------------------- */
499
500 /* Note: enough space needs to be allocated here so that routines called by
501 * cholmod_analyze do not reallocate the space.
502 */
503
504 /* s = 6*n + uncol */
505 s = CHOLMOD(mult_size_t) (n, 6, &ok) ;
506 s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
507 if (!ok)
508 {
509 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
510 return (NULL) ;
511 }
512
513 CHOLMOD(allocate_work) (n, s, 0, Common) ;
514 if (Common->status < CHOLMOD_OK)
515 {
516 return (NULL) ; /* out of memory */
517 }
518 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
519
520 /* ensure that subsequent routines, called by cholmod_analyze, do not
521 * reallocate any workspace. This is set back to FALSE in the
522 * FREE_WORKSPACE_AND_RETURN macro, which is the only way this function
523 * returns to its caller. */
524 Common->no_workspace_reallocate = TRUE ;
525
526 /* Use the last 4*n Int's in Iwork for Parent, First, Level, and Post, since
527 * other CHOLMOD routines will use the first 2n+uncol space. The ordering
528 * routines (cholmod_amd, cholmod_colamd, cholmod_ccolamd, cholmod_metis)
529 * are an exception. They can use all 6n + ncol space, since the contents
530 * of Parent, First, Level, and Post are not needed across calls to those
531 * routines. */
532 Work4n = Common->Iwork ;
533 Work4n += 2*((size_t) n) + uncol ;
534 Parent = Work4n ;
535 First = Work4n + n ;
536 Level = Work4n + 2*((size_t) n) ;
537 Post = Work4n + 3*((size_t) n) ;
538
539 /* note that this assignment means that cholmod_nested_dissection,
540 * cholmod_ccolamd, and cholmod_camd can use only the first 4n+uncol
541 * space in Common->Iwork */
542 Cmember = Post ;
543 CParent = Level ;
544
545 /* ---------------------------------------------------------------------- */
546 /* allocate more workspace, and an empty simplicial symbolic factor */
547 /* ---------------------------------------------------------------------- */
548
549 L = CHOLMOD(allocate_factor) (n, Common) ;
550 Lparent = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
551 Perm = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
552 ColCount = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
553 if (Common->status < CHOLMOD_OK)
554 {
555 /* out of memory */
556 FREE_WORKSPACE_AND_RETURN ;
557 }
558 Lperm = L->Perm ;
559 Lcolcount = L->ColCount ;
560 Common->anz = EMPTY ;
561
562 /* ---------------------------------------------------------------------- */
563 /* try all the requested ordering options and backup to AMD if needed */
564 /* ---------------------------------------------------------------------- */
565
566 /* turn off error handling [ */
567 Common->try_catch = TRUE ;
568
569 for (method = 0 ; method <= nmethods ; method++)
570 {
571
572 /* ------------------------------------------------------------------ */
573 /* determine the method to try */
574 /* ------------------------------------------------------------------ */
575
576 Common->fl = EMPTY ;
577 Common->lnz = EMPTY ;
578 skip_analysis = FALSE ;
579
580 if (method == nmethods)
581 {
582 /* All methods failed: backup to AMD */
583 if (Common->selected == EMPTY && amd_backup)
584 {
585 PRINT1 (("All methods requested failed: backup to AMD\n")) ;
586 ordering = CHOLMOD_AMD ;
587 }
588 else
589 {
590 break ;
591 }
592 }
593 else
594 {
595 ordering = Common->method [method].ordering ;
596 }
597 Common->current = method ;
598 PRINT1 (("method "ID": Try method: "ID"\n", method, ordering)) ;
599
600 /* ------------------------------------------------------------------ */
601 /* find the fill-reducing permutation */
602 /* ------------------------------------------------------------------ */
603
604 if (ordering == CHOLMOD_NATURAL)
605 {
606
607 /* -------------------------------------------------------------- */
608 /* natural ordering */
609 /* -------------------------------------------------------------- */
610
611 for (k = 0 ; k < n ; k++)
612 {
613 Perm [k] = k ;
614 }
615
616 }
617 else if (ordering == CHOLMOD_GIVEN)
618 {
619
620 /* -------------------------------------------------------------- */
621 /* use given ordering of A, if provided */
622 /* -------------------------------------------------------------- */
623
624 if (UserPerm == NULL)
625 {
626 /* this is not an error condition */
627 PRINT1 (("skip, no user perm given\n")) ;
628 continue ;
629 }
630 for (k = 0 ; k < n ; k++)
631 {
632 /* UserPerm is checked in cholmod_ptranspose */
633 Perm [k] = UserPerm [k] ;
634 }
635
636 }
637 else if (ordering == CHOLMOD_AMD)
638 {
639
640 /* -------------------------------------------------------------- */
641 /* AMD ordering of A, A*A', or A(:,f)*A(:,f)' */
642 /* -------------------------------------------------------------- */
643
644 amd_backup = FALSE ; /* no need to try AMD twice ... */
645 CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
646 skip_analysis = TRUE ;
647
648 }
649 else if (ordering == CHOLMOD_COLAMD)
650 {
651
652 /* -------------------------------------------------------------- */
653 /* AMD for symmetric case, COLAMD for A*A' or A(:,f)*A(:,f)' */
654 /* -------------------------------------------------------------- */
655
656 if (A->stype)
657 {
658 CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
659 skip_analysis = TRUE ;
660 }
661 else
662 {
663 /* Alternative:
664 CHOLMOD(ccolamd) (A, fset, fsize, NULL, Perm, Common) ;
665 */
666 /* do not postorder, it is done later, below */
667 /* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1)*/
668 CHOLMOD(colamd) (A, fset, fsize, FALSE, Perm, Common) ;
669 }
670
671 }
672 else if (ordering == CHOLMOD_METIS)
673 {
674
675 /* -------------------------------------------------------------- */
676 /* use METIS_NodeND directly (via a CHOLMOD wrapper) */
677 /* -------------------------------------------------------------- */
678
679 #ifndef NPARTITION
680 /* postorder parameter is false, because it will be later, below */
681 /* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1) */
682 Common->called_nd = TRUE ;
683 CHOLMOD(metis) (A, fset, fsize, FALSE, Perm, Common) ;
684 #else
685 Common->status = CHOLMOD_NOT_INSTALLED ;
686 #endif
687
688 }
689 else if (ordering == CHOLMOD_NESDIS)
690 {
691
692 /* -------------------------------------------------------------- */
693 /* use CHOLMOD's nested dissection */
694 /* -------------------------------------------------------------- */
695
696 /* this method is based on METIS' node bissection routine
697 * (METIS_ComputeVertexSeparator). In contrast to METIS_NodeND,
698 * it calls CAMD or CCOLAMD on the whole graph, instead of MMD
699 * on just the leaves. */
700 #ifndef NPARTITION
701 /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow) */
702 Common->called_nd = TRUE ;
703 CHOLMOD(nested_dissection) (A, fset, fsize, Perm, CParent, Cmember,
704 Common) ;
705 #else
706 Common->status = CHOLMOD_NOT_INSTALLED ;
707 #endif
708
709 }
710 else
711 {
712
713 /* -------------------------------------------------------------- */
714 /* invalid ordering method */
715 /* -------------------------------------------------------------- */
716
717 Common->status = CHOLMOD_INVALID ;
718 PRINT1 (("No such ordering: "ID"\n", ordering)) ;
719 }
720
721 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
722
723 if (Common->status < CHOLMOD_OK)
724 {
725 /* out of memory, or method failed */
726 status = MIN (status, Common->status) ;
727 Common->status = CHOLMOD_OK ;
728 continue ;
729 }
730
731 /* ------------------------------------------------------------------ */
732 /* analyze the ordering */
733 /* ------------------------------------------------------------------ */
734
735 if (!skip_analysis)
736 {
737 if (!CHOLMOD(analyze_ordering) (A, ordering, Perm, fset, fsize,
738 Parent, Post, ColCount, First, Level, Common))
739 {
740 /* ordering method failed; clear status and try next method */
741 status = MIN (status, Common->status) ;
742 Common->status = CHOLMOD_OK ;
743 continue ;
744 }
745 }
746
747 ASSERT (Common->fl >= 0 && Common->lnz >= 0) ;
748 Common->method [method].fl = Common->fl ;
749 Common->method [method].lnz = Common->lnz ;
750 PRINT1 (("lnz %g fl %g\n", Common->lnz, Common->fl)) ;
751
752 /* ------------------------------------------------------------------ */
753 /* pick the best method */
754 /* ------------------------------------------------------------------ */
755
756 /* fl.pt. compare, but lnz can never be NaN */
757 if (Common->selected == EMPTY || Common->lnz < lnz_best)
758 {
759 Common->selected = method ;
760 PRINT1 (("this is best so far, method "ID"\n", method)) ;
761 L->ordering = ordering ;
762 lnz_best = Common->lnz ;
763 for (k = 0 ; k < n ; k++)
764 {
765 Lperm [k] = Perm [k] ;
766 }
767 /* save the results of cholmod_analyze_ordering, if it was called */
768 skip_best = skip_analysis ;
769 if (!skip_analysis)
770 {
771 /* save the column count; becomes permanent part of L */
772 for (k = 0 ; k < n ; k++)
773 {
774 Lcolcount [k] = ColCount [k] ;
775 }
776 /* Parent is needed for weighted postordering and for supernodal
777 * analysis. Does not become a permanent part of L */
778 for (k = 0 ; k < n ; k++)
779 {
780 Lparent [k] = Parent [k] ;
781 }
782 }
783 }
784
785 /* ------------------------------------------------------------------ */
786 /* determine if METIS is to be skipped */
787 /* ------------------------------------------------------------------ */
788
789 if (default_strategy && ordering == CHOLMOD_AMD)
790 {
791 if ((Common->fl < 500 * Common->lnz) ||
792 (Common->lnz < 5 * Common->anz))
793 {
794 /* AMD found an ordering with less than 500 flops per nonzero in
795 * L, or one with a fill-in ratio (nnz(L)/nnz(A)) of less than
796 * 5. This is pretty good, and it's unlikely that METIS will do
797 * better (this heuristic is based on tests on all symmetric
798 * positive definite matrices in the UF sparse matrix
799 * collection, and it works well across a wide range of
800 * problems). METIS can take much more time than AMD. */
801 break ;
802 }
803 }
804 }
805
806 /* turn error printing back on ] */
807 Common->try_catch = FALSE ;
808
809 /* ---------------------------------------------------------------------- */
810 /* return if no ordering method succeeded */
811 /* ---------------------------------------------------------------------- */
812
813 if (Common->selected == EMPTY)
814 {
815 /* All methods failed.
816 * If two or more methods failed, they may have failed for different
817 * reasons. Both would clear Common->status and skip to the next
818 * method. Common->status needs to be restored here to the worst error
819 * obtained in any of the methods. CHOLMOD_INVALID is worse
820 * than CHOLMOD_OUT_OF_MEMORY, since the former implies something may
821 * be wrong with the user's input. CHOLMOD_OUT_OF_MEMORY is simply an
822 * indication of lack of resources. */
823 if (status >= CHOLMOD_OK)
824 {
825 /* this can occur if nmethods = 1, ordering = CHOLMOD_GIVEN,
826 but UserPerm is NULL */
827 status = CHOLMOD_INVALID ;
828 }
829 ERROR (status, "all methods failed") ;
830 FREE_WORKSPACE_AND_RETURN ;
831 }
832
833 /* ---------------------------------------------------------------------- */
834 /* do the analysis for AMD, if skipped */
835 /* ---------------------------------------------------------------------- */
836
837 Common->fl = Common->method [Common->selected].fl ;
838 Common->lnz = Common->method [Common->selected].lnz ;
839 ASSERT (Common->lnz >= 0) ;
840
841 if (skip_best)
842 {
843 if (!CHOLMOD(analyze_ordering) (A, L->ordering, Lperm, fset, fsize,
844 Lparent, Post, Lcolcount, First, Level, Common))
845 {
846 /* out of memory, or method failed */
847 FREE_WORKSPACE_AND_RETURN ;
848 }
849 }
850
851 /* ---------------------------------------------------------------------- */
852 /* postorder the etree, weighted by the column counts */
853 /* ---------------------------------------------------------------------- */
854
855 if (Common->postorder)
856 {
857 /* combine the fill-reducing ordering with the weighted postorder */
858 /* workspace: Iwork (2*nrow) */
859 if (CHOLMOD(postorder) (Lparent, n, Lcolcount, Post, Common) == n)
860 {
861 /* use First and Level as workspace [ */
862 Int *Wi = First, *InvPost = Level ;
863 Int newchild, oldchild, newparent, oldparent ;
864
865 for (k = 0 ; k < n ; k++)
866 {
867 Wi [k] = Lperm [Post [k]] ;
868 }
869 for (k = 0 ; k < n ; k++)
870 {
871 Lperm [k] = Wi [k] ;
872 }
873
874 for (k = 0 ; k < n ; k++)
875 {
876 Wi [k] = Lcolcount [Post [k]] ;
877 }
878 for (k = 0 ; k < n ; k++)
879 {
880 Lcolcount [k] = Wi [k] ;
881 }
882 for (k = 0 ; k < n ; k++)
883 {
884 InvPost [Post [k]] = k ;
885 }
886
887 /* updated Lparent needed only for supernodal case */
888 for (newchild = 0 ; newchild < n ; newchild++)
889 {
890 oldchild = Post [newchild] ;
891 oldparent = Lparent [oldchild] ;
892 newparent = (oldparent == EMPTY) ? EMPTY : InvPost [oldparent] ;
893 Wi [newchild] = newparent ;
894 }
895 for (k = 0 ; k < n ; k++)
896 {
897 Lparent [k] = Wi [k] ;
898 }
899 /* done using Iwork as workspace ] */
900
901 /* L is now postordered, no longer in natural ordering */
902 if (L->ordering == CHOLMOD_NATURAL)
903 {
904 L->ordering = CHOLMOD_POSTORDERED ;
905 }
906 }
907 }
908
909 /* ---------------------------------------------------------------------- */
910 /* supernodal analysis, if requested or if selected automatically */
911 /* ---------------------------------------------------------------------- */
912
913 #ifndef NSUPERNODAL
914 if (Common->supernodal > CHOLMOD_AUTO
915 || (Common->supernodal == CHOLMOD_AUTO &&
916 Common->lnz > 0 &&
917 (Common->fl / Common->lnz) >= Common->supernodal_switch))
918 {
919 cholmod_sparse *S, *F, *A2, *A1 ;
920
921 permute_matrices (A, L->ordering, Lperm, fset, fsize, TRUE,
922 &A1, &A2, &S, &F, Common) ;
923
924 /* workspace: Flag (nrow), Head (nrow), Iwork (5*nrow) */
925 CHOLMOD(super_symbolic2) (for_whom, S, F, Lparent, L, Common) ;
926 PRINT1 (("status %d\n", Common->status)) ;
927
928 CHOLMOD(free_sparse) (&A1, Common) ;
929 CHOLMOD(free_sparse) (&A2, Common) ;
930 }
931 #endif
932
933 /* ---------------------------------------------------------------------- */
934 /* free temporary matrices and workspace, and return result L */
935 /* ---------------------------------------------------------------------- */
936
937 FREE_WORKSPACE_AND_RETURN ;
938 }
939 #endif
940