1 /* ========================================================================== */
2 /* === Cholesky/cholmod_etree =============================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Compute the elimination tree of A or A'*A
10 *
11 * In the symmetric case, the upper triangular part of A is used. Entries not
12 * in this part of the matrix are ignored. Computing the etree of a symmetric
13 * matrix from just its lower triangular entries is not supported.
14 *
15 * In the unsymmetric case, all of A is used, and the etree of A'*A is computed.
16 *
17 * References:
18 *
19 * J. Liu, "A compact row storage scheme for Cholesky factors", ACM Trans.
20 * Math. Software, vol 12, 1986, pp. 127-148.
21 *
22 * J. Liu, "The role of elimination trees in sparse factorization", SIAM J.
23 * Matrix Analysis & Applic., vol 11, 1990, pp. 134-172.
24 *
25 * J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
26 * sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
27 *
28 * workspace: symmetric: Iwork (nrow), unsymmetric: Iwork (nrow+ncol)
29 *
30 * Supports any xtype (pattern, real, complex, or zomplex)
31 */
32
33 #ifndef NCHOLESKY
34
35 #include "cholmod_internal.h"
36 #include "cholmod_cholesky.h"
37
38 /* ========================================================================== */
39 /* === update_etree ========================================================= */
40 /* ========================================================================== */
41
update_etree(Int k,Int i,Int Parent[],Int Ancestor[])42 static void update_etree
43 (
44 /* inputs, not modified */
45 Int k, /* process the edge (k,i) in the input graph */
46 Int i,
47 /* inputs, modified on output */
48 Int Parent [ ], /* Parent [t] = p if p is the parent of t */
49 Int Ancestor [ ] /* Ancestor [t] is the ancestor of node t in the
50 partially-constructed etree */
51 )
52 {
53 Int a ;
54 for ( ; ; ) /* traverse the path from k to the root of the tree */
55 {
56 a = Ancestor [k] ;
57 if (a == i)
58 {
59 /* final ancestor reached; no change to tree */
60 return ;
61 }
62 /* perform path compression */
63 Ancestor [k] = i ;
64 if (a == EMPTY)
65 {
66 /* final ancestor undefined; this is a new edge in the tree */
67 Parent [k] = i ;
68 return ;
69 }
70 /* traverse up to the ancestor of k */
71 k = a ;
72 }
73 }
74
75 /* ========================================================================== */
76 /* === cholmod_etree ======================================================== */
77 /* ========================================================================== */
78
79 /* Find the elimination tree of A or A'*A */
80
CHOLMOD(etree)81 int CHOLMOD(etree)
82 (
83 /* ---- input ---- */
84 cholmod_sparse *A,
85 /* ---- output --- */
86 Int *Parent, /* size ncol. Parent [j] = p if p is the parent of j */
87 /* --------------- */
88 cholmod_common *Common
89 )
90 {
91 Int *Ap, *Ai, *Anz, *Ancestor, *Prev, *Iwork ;
92 Int i, j, jprev, p, pend, nrow, ncol, packed, stype ;
93 size_t s ;
94 int ok = TRUE ;
95
96 /* ---------------------------------------------------------------------- */
97 /* check inputs */
98 /* ---------------------------------------------------------------------- */
99
100 RETURN_IF_NULL_COMMON (FALSE) ;
101 RETURN_IF_NULL (A, FALSE) ;
102 RETURN_IF_NULL (Parent, FALSE) ;
103 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
104 Common->status = CHOLMOD_OK ;
105
106 /* ---------------------------------------------------------------------- */
107 /* allocate workspace */
108 /* ---------------------------------------------------------------------- */
109
110 stype = A->stype ;
111
112 /* s = A->nrow + (stype ? 0 : A->ncol) */
113 s = CHOLMOD(add_size_t) (A->nrow, (stype ? 0 : A->ncol), &ok) ;
114 if (!ok)
115 {
116 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
117 return (FALSE) ;
118 }
119
120 CHOLMOD(allocate_work) (0, s, 0, Common) ;
121 if (Common->status < CHOLMOD_OK)
122 {
123 return (FALSE) ; /* out of memory */
124 }
125
126 ASSERT (CHOLMOD(dump_sparse) (A, "etree", Common) >= 0) ;
127 Iwork = Common->Iwork ;
128
129 /* ---------------------------------------------------------------------- */
130 /* get inputs */
131 /* ---------------------------------------------------------------------- */
132
133 ncol = A->ncol ; /* the number of columns of A */
134 nrow = A->nrow ; /* the number of rows of A */
135 Ap = A->p ; /* size ncol+1, column pointers for A */
136 Ai = A->i ; /* the row indices of A */
137 Anz = A->nz ; /* number of nonzeros in each column of A */
138 packed = A->packed ;
139 Ancestor = Iwork ; /* size ncol (i/i/l) */
140
141 for (j = 0 ; j < ncol ; j++)
142 {
143 Parent [j] = EMPTY ;
144 Ancestor [j] = EMPTY ;
145 }
146
147 /* ---------------------------------------------------------------------- */
148 /* compute the etree */
149 /* ---------------------------------------------------------------------- */
150
151 if (stype > 0)
152 {
153
154 /* ------------------------------------------------------------------ */
155 /* symmetric (upper) case: compute etree (A) */
156 /* ------------------------------------------------------------------ */
157
158 for (j = 0 ; j < ncol ; j++)
159 {
160 /* for each row i in column j of triu(A), excluding the diagonal */
161 p = Ap [j] ;
162 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
163 for ( ; p < pend ; p++)
164 {
165 i = Ai [p] ;
166 if (i < j)
167 {
168 update_etree (i, j, Parent, Ancestor) ;
169 }
170 }
171 }
172
173 }
174 else if (stype == 0)
175 {
176
177 /* ------------------------------------------------------------------ */
178 /* unsymmetric case: compute etree (A'*A) */
179 /* ------------------------------------------------------------------ */
180
181 Prev = Iwork + ncol ; /* size nrow (i/i/l) */
182 for (i = 0 ; i < nrow ; i++)
183 {
184 Prev [i] = EMPTY ;
185 }
186 for (j = 0 ; j < ncol ; j++)
187 {
188 /* for each row i in column j of A */
189 p = Ap [j] ;
190 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
191 for ( ; p < pend ; p++)
192 {
193 /* a graph is constructed dynamically with one path per row
194 * of A. If the ith row of A contains column indices
195 * (j1,j2,j3,j4) then the new graph has edges (j1,j2), (j2,j3),
196 * and (j3,j4). When at node i of this path-graph, all edges
197 * (jprev,j) are considered, where jprev<j */
198 i = Ai [p] ;
199 jprev = Prev [i] ;
200 if (jprev != EMPTY)
201 {
202 update_etree (jprev, j, Parent, Ancestor) ;
203 }
204 Prev [i] = j ;
205 }
206 }
207
208 }
209 else
210 {
211
212 /* ------------------------------------------------------------------ */
213 /* symmetric case with lower triangular part not supported */
214 /* ------------------------------------------------------------------ */
215
216 ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
217 return (FALSE) ;
218 }
219
220 ASSERT (CHOLMOD(dump_parent) (Parent, ncol, "Parent", Common)) ;
221 return (TRUE) ;
222 }
223 #endif
224