1 /* ========================================================================== */
2 /* === Cholesky/cholmod_colamd ============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* CHOLMOD interface to the COLAMD ordering routine (version 2.4 or later).
10 * Finds a permutation p such that the Cholesky factorization of PAA'P' is
11 * sparser than AA' using colamd. If the postorder input parameter is TRUE,
12 * the column etree is found and postordered, and the colamd ordering is then
13 * combined with its postordering. A must be unsymmetric.
14 *
15 * There can be no duplicate entries in f.
16 * f can be length 0 to n if A is m-by-n.
17 *
18 * workspace: Iwork (4*nrow+ncol), Head (nrow+1), Flag (nrow)
19 * Allocates a copy of its input matrix, which
20 * is then used as CCOLAMD's workspace.
21 *
22 * Supports any xtype (pattern, real, complex, or zomplex)
23 */
24
25 #ifndef NCHOLESKY
26
27 #include "cholmod_internal.h"
28 #include "colamd.h"
29 #include "cholmod_cholesky.h"
30
31 #if (!defined (COLAMD_VERSION) || (COLAMD_VERSION < COLAMD_VERSION_CODE (2,5)))
32 #error "COLAMD v2.5 or later is required"
33 #endif
34
35 /* ========================================================================== */
36 /* === cholmod_colamd ======================================================= */
37 /* ========================================================================== */
38
CHOLMOD(colamd)39 int CHOLMOD(colamd)
40 (
41 /* ---- input ---- */
42 cholmod_sparse *A, /* matrix to order */
43 Int *fset, /* subset of 0:(A->ncol)-1 */
44 size_t fsize, /* size of fset */
45 int postorder, /* if TRUE, follow with a coletree postorder */
46 /* ---- output --- */
47 Int *Perm, /* size A->nrow, output permutation */
48 /* --------------- */
49 cholmod_common *Common
50 )
51 {
52 double knobs [COLAMD_KNOBS] ;
53 cholmod_sparse *C ;
54 Int *NewPerm, *Parent, *Post, *Work2n ;
55 Int k, nrow, ncol ;
56 size_t s, alen ;
57 int ok = TRUE ;
58
59 /* ---------------------------------------------------------------------- */
60 /* check inputs */
61 /* ---------------------------------------------------------------------- */
62
63 RETURN_IF_NULL_COMMON (FALSE) ;
64 RETURN_IF_NULL (A, FALSE) ;
65 RETURN_IF_NULL (Perm, FALSE) ;
66 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
67 if (A->stype != 0)
68 {
69 ERROR (CHOLMOD_INVALID, "matrix must be unsymmetric") ;
70 return (FALSE) ;
71 }
72 Common->status = CHOLMOD_OK ;
73
74 /* ---------------------------------------------------------------------- */
75 /* get inputs */
76 /* ---------------------------------------------------------------------- */
77
78 nrow = A->nrow ;
79 ncol = A->ncol ;
80
81 /* ---------------------------------------------------------------------- */
82 /* allocate workspace */
83 /* ---------------------------------------------------------------------- */
84
85 /* Note: this is less than the space used in cholmod_analyze, so if
86 * cholmod_colamd is being called by that routine, no space will be
87 * allocated.
88 */
89
90 /* s = 4*nrow + ncol */
91 s = CHOLMOD(mult_size_t) (nrow, 4, &ok) ;
92 s = CHOLMOD(add_size_t) (s, ncol, &ok) ;
93
94 #ifdef LONG
95 alen = colamd_l_recommended (A->nzmax, ncol, nrow) ;
96 colamd_l_set_defaults (knobs) ;
97 #else
98 alen = colamd_recommended (A->nzmax, ncol, nrow) ;
99 colamd_set_defaults (knobs) ;
100 #endif
101
102 if (!ok || alen == 0)
103 {
104 ERROR (CHOLMOD_TOO_LARGE, "matrix invalid or too large") ;
105 return (FALSE) ;
106 }
107
108 CHOLMOD(allocate_work) (0, s, 0, Common) ;
109 if (Common->status < CHOLMOD_OK)
110 {
111 return (FALSE) ;
112 }
113
114 /* ---------------------------------------------------------------------- */
115 /* allocate COLAMD workspace */
116 /* ---------------------------------------------------------------------- */
117
118 C = CHOLMOD(allocate_sparse) (ncol, nrow, alen, TRUE, TRUE, 0,
119 CHOLMOD_PATTERN, Common) ;
120
121 /* ---------------------------------------------------------------------- */
122 /* copy (and transpose) the input matrix A into the colamd workspace */
123 /* ---------------------------------------------------------------------- */
124
125 /* C = A (:,f)', which also packs A if needed. */
126 /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset) */
127 ok = CHOLMOD(transpose_unsym) (A, 0, NULL, fset, fsize, C, Common) ;
128
129 /* ---------------------------------------------------------------------- */
130 /* order the matrix (destroys the contents of C->i and C->p) */
131 /* ---------------------------------------------------------------------- */
132
133 /* get parameters */
134 if (Common->current < 0 || Common->current >= CHOLMOD_MAXMETHODS)
135 {
136 /* this is the CHOLMOD default, not the COLAMD default */
137 knobs [COLAMD_DENSE_ROW] = -1 ;
138 }
139 else
140 {
141 /* get the knobs from the Common parameters */
142 knobs [COLAMD_DENSE_COL] = Common->method[Common->current].prune_dense ;
143 knobs [COLAMD_DENSE_ROW] = Common->method[Common->current].prune_dense2;
144 knobs [COLAMD_AGGRESSIVE] = Common->method[Common->current].aggressive ;
145 }
146
147 if (ok)
148 {
149 Int *Cp ;
150 Int stats [COLAMD_STATS] ;
151 Cp = C->p ;
152
153 #ifdef LONG
154 colamd_l (ncol, nrow, alen, C->i, Cp, knobs, stats) ;
155 #else
156 colamd (ncol, nrow, alen, C->i, Cp, knobs, stats) ;
157 #endif
158
159 ok = stats [COLAMD_STATUS] ;
160 ok = (ok == COLAMD_OK || ok == COLAMD_OK_BUT_JUMBLED) ;
161 /* permutation returned in C->p, if the ordering succeeded */
162 for (k = 0 ; k < nrow ; k++)
163 {
164 Perm [k] = Cp [k] ;
165 }
166 }
167
168 CHOLMOD(free_sparse) (&C, Common) ;
169
170 /* ---------------------------------------------------------------------- */
171 /* column etree postordering */
172 /* ---------------------------------------------------------------------- */
173
174 if (postorder)
175 {
176 /* use the last 2*n space in Iwork for Parent and Post */
177 Work2n = Common->Iwork ;
178 Work2n += 2*((size_t) nrow) + ncol ;
179 Parent = Work2n ; /* size nrow (i/i/l) */
180 Post = Work2n + nrow ; /* size nrow (i/i/l) */
181
182 /* workspace: Iwork (2*nrow+ncol), Flag (nrow), Head (nrow+1) */
183 ok = ok && CHOLMOD(analyze_ordering) (A, CHOLMOD_COLAMD, Perm, fset,
184 fsize, Parent, Post, NULL, NULL, NULL, Common) ;
185
186 /* combine the colamd permutation with its postordering */
187 if (ok)
188 {
189 NewPerm = Common->Iwork ; /* size nrow (i/i/l) */
190 for (k = 0 ; k < nrow ; k++)
191 {
192 NewPerm [k] = Perm [Post [k]] ;
193 }
194 for (k = 0 ; k < nrow ; k++)
195 {
196 Perm [k] = NewPerm [k] ;
197 }
198 }
199 }
200
201 return (ok) ;
202 }
203 #endif
204