|
0
|
1 |
|
|
2 /* Elimination tree computation and layout routines */ |
|
|
3 |
|
|
4 #include <stdio.h> |
|
|
5 #include <malloc.h> |
|
|
6 #include <stdlib.h> |
|
|
7 #include "util.h" |
|
|
8 |
|
|
9 /* |
|
|
10 * Implementation of disjoint set union routines. |
|
|
11 * Elements are integers in 0..n-1, and the |
|
|
12 * names of the sets themselves are of type int. |
|
|
13 * |
|
|
14 * Calls are: |
|
|
15 * initialize_disjoint_sets (n) initial call. |
|
|
16 * s = make_set (i) returns a set containing only i. |
|
|
17 * s = link (t, u) returns s = t union u, destroying t and u. |
|
|
18 * s = find (i) return name of set containing i. |
|
|
19 * finalize_disjoint_sets final call. |
|
|
20 * |
|
|
21 * This implementation uses path compression but not weighted union. |
|
|
22 * See Tarjan's book for details. |
|
|
23 * John Gilbert, CMI, 1987. |
|
|
24 * |
|
|
25 * Implemented path-halving by XL 7/5/95. |
|
|
26 */ |
|
|
27 |
|
|
28 static int *pp; /* parent array for sets */ |
|
|
29 |
|
|
30 static |
|
|
31 int *mxCallocInt(int n) |
|
|
32 { |
|
|
33 register int i; |
|
|
34 int *buf; |
|
|
35 |
|
|
36 buf = (int *) SUPERLU_MALLOC( n * sizeof(int) ); |
|
|
37 if ( !buf ) { |
|
|
38 ABORT("SUPERLU_MALLOC fails for buf in mxCallocInt()"); |
|
|
39 } |
|
|
40 for (i = 0; i < n; i++) buf[i] = 0; |
|
|
41 return (buf); |
|
|
42 } |
|
|
43 |
|
|
44 static |
|
|
45 void initialize_disjoint_sets ( |
|
|
46 int n |
|
|
47 ) |
|
|
48 { |
|
|
49 pp = mxCallocInt(n); |
|
|
50 } |
|
|
51 |
|
|
52 |
|
|
53 static |
|
|
54 int make_set ( |
|
|
55 int i |
|
|
56 ) |
|
|
57 { |
|
|
58 pp[i] = i; |
|
|
59 return i; |
|
|
60 } |
|
|
61 |
|
|
62 |
|
|
63 static |
|
|
64 int link ( |
|
|
65 int s, |
|
|
66 int t |
|
|
67 ) |
|
|
68 { |
|
|
69 pp[s] = t; |
|
|
70 return t; |
|
|
71 } |
|
|
72 |
|
|
73 |
|
|
74 /* PATH HALVING */ |
|
|
75 static |
|
|
76 int find (int i) |
|
|
77 { |
|
|
78 register int p, gp; |
|
|
79 |
|
|
80 p = pp[i]; |
|
|
81 gp = pp[p]; |
|
|
82 while (gp != p) { |
|
|
83 pp[i] = gp; |
|
|
84 i = gp; |
|
|
85 p = pp[i]; |
|
|
86 gp = pp[p]; |
|
|
87 } |
|
|
88 return (p); |
|
|
89 } |
|
|
90 |
|
|
91 #if 0 |
|
|
92 /* PATH COMPRESSION */ |
|
|
93 static |
|
|
94 int find ( |
|
|
95 int i |
|
|
96 ) |
|
|
97 { |
|
|
98 if (pp[i] != i) |
|
|
99 pp[i] = find (pp[i]); |
|
|
100 return pp[i]; |
|
|
101 } |
|
|
102 #endif |
|
|
103 |
|
|
104 static |
|
|
105 void finalize_disjoint_sets ( |
|
|
106 void |
|
|
107 ) |
|
|
108 { |
|
|
109 SUPERLU_FREE(pp); |
|
|
110 } |
|
|
111 |
|
|
112 |
|
|
113 /* |
|
|
114 * Find the elimination tree for A'*A. |
|
|
115 * This uses something similar to Liu's algorithm. |
|
|
116 * It runs in time O(nz(A)*log n) and does not form A'*A. |
|
|
117 * |
|
|
118 * Input: |
|
|
119 * Sparse matrix A. Numeric values are ignored, so any |
|
|
120 * explicit zeros are treated as nonzero. |
|
|
121 * Output: |
|
|
122 * Integer array of parents representing the elimination |
|
|
123 * tree of the symbolic product A'*A. Each vertex is a |
|
|
124 * column of A, and nc means a root of the elimination forest. |
|
|
125 * |
|
|
126 * John R. Gilbert, Xerox, 10 Dec 1990 |
|
|
127 * Based on code by JRG dated 1987, 1988, and 1990. |
|
|
128 */ |
|
|
129 |
|
|
130 /* |
|
|
131 * Nonsymmetric elimination tree |
|
|
132 */ |
|
|
133 int |
|
|
134 sp_coletree( |
|
|
135 int *acolst, int *acolend, /* column start and end past 1 */ |
|
|
136 int *arow, /* row indices of A */ |
|
|
137 int nr, int nc, /* dimension of A */ |
|
|
138 int *parent /* parent in elim tree */ |
|
|
139 ) |
|
|
140 { |
|
|
141 int *root; /* root of subtee of etree */ |
|
|
142 int *firstcol; /* first nonzero col in each row*/ |
|
|
143 int rset, cset; |
|
|
144 int row, col; |
|
|
145 int rroot; |
|
|
146 int p; |
|
|
147 |
|
|
148 root = mxCallocInt (nc); |
|
|
149 initialize_disjoint_sets (nc); |
|
|
150 |
|
|
151 /* Compute firstcol[row] = first nonzero column in row */ |
|
|
152 |
|
|
153 firstcol = mxCallocInt (nr); |
|
|
154 for (row = 0; row < nr; firstcol[row++] = nc); |
|
|
155 for (col = 0; col < nc; col++) |
|
|
156 for (p = acolst[col]; p < acolend[col]; p++) { |
|
|
157 row = arow[p]; |
|
|
158 firstcol[row] = MIN(firstcol[row], col); |
|
|
159 } |
|
|
160 |
|
|
161 /* Compute etree by Liu's algorithm for symmetric matrices, |
|
|
162 except use (firstcol[r],c) in place of an edge (r,c) of A. |
|
|
163 Thus each row clique in A'*A is replaced by a star |
|
|
164 centered at its first vertex, which has the same fill. */ |
|
|
165 |
|
|
166 for (col = 0; col < nc; col++) { |
|
|
167 cset = make_set (col); |
|
|
168 root[cset] = col; |
|
|
169 parent[col] = nc; /* Matlab */ |
|
|
170 for (p = acolst[col]; p < acolend[col]; p++) { |
|
|
171 row = firstcol[arow[p]]; |
|
|
172 if (row >= col) continue; |
|
|
173 rset = find (row); |
|
|
174 rroot = root[rset]; |
|
|
175 if (rroot != col) { |
|
|
176 parent[rroot] = col; |
|
|
177 cset = link (cset, rset); |
|
|
178 root[cset] = col; |
|
|
179 } |
|
|
180 } |
|
|
181 } |
|
|
182 |
|
|
183 SUPERLU_FREE (root); |
|
|
184 SUPERLU_FREE (firstcol); |
|
|
185 finalize_disjoint_sets (); |
|
|
186 return 0; |
|
|
187 } |
|
|
188 |
|
|
189 /* |
|
|
190 * q = TreePostorder (n, p); |
|
|
191 * |
|
|
192 * Postorder a tree. |
|
|
193 * Input: |
|
|
194 * p is a vector of parent pointers for a forest whose |
|
|
195 * vertices are the integers 0 to n-1; p[root]==n. |
|
|
196 * Output: |
|
|
197 * q is a vector indexed by 0..n-1 such that q[i] is the |
|
|
198 * i-th vertex in a postorder numbering of the tree. |
|
|
199 * |
|
|
200 * ( 2/7/95 modified by X.Li: |
|
|
201 * q is a vector indexed by 0:n-1 such that vertex i is the |
|
|
202 * q[i]-th vertex in a postorder numbering of the tree. |
|
|
203 * That is, this is the inverse of the previous q. ) |
|
|
204 * |
|
|
205 * In the child structure, lower-numbered children are represented |
|
|
206 * first, so that a tree which is already numbered in postorder |
|
|
207 * will not have its order changed. |
|
|
208 * |
|
|
209 * Written by John Gilbert, Xerox, 10 Dec 1990. |
|
|
210 * Based on code written by John Gilbert at CMI in 1987. |
|
|
211 */ |
|
|
212 |
|
|
213 static int *first_kid, *next_kid; /* Linked list of children. */ |
|
|
214 static int *post, postnum; |
|
|
215 |
|
|
216 static |
|
|
217 /* |
|
|
218 * Depth-first search from vertex v. |
|
|
219 */ |
|
|
220 void etdfs ( |
|
|
221 int v |
|
|
222 ) |
|
|
223 { |
|
|
224 int w; |
|
|
225 |
|
|
226 for (w = first_kid[v]; w != -1; w = next_kid[w]) { |
|
|
227 etdfs (w); |
|
|
228 } |
|
|
229 /* post[postnum++] = v; in Matlab */ |
|
|
230 post[v] = postnum++; /* Modified by X.Li on 2/14/95 */ |
|
|
231 } |
|
|
232 |
|
|
233 |
|
|
234 /* |
|
|
235 * Post order a tree |
|
|
236 */ |
|
|
237 int *TreePostorder( |
|
|
238 int n, |
|
|
239 int *parent |
|
|
240 ) |
|
|
241 { |
|
|
242 int v, dad; |
|
|
243 |
|
|
244 /* Allocate storage for working arrays and results */ |
|
|
245 first_kid = mxCallocInt (n+1); |
|
|
246 next_kid = mxCallocInt (n+1); |
|
|
247 post = mxCallocInt (n+1); |
|
|
248 |
|
|
249 /* Set up structure describing children */ |
|
|
250 for (v = 0; v <= n; first_kid[v++] = -1); |
|
|
251 for (v = n-1; v >= 0; v--) { |
|
|
252 dad = parent[v]; |
|
|
253 next_kid[v] = first_kid[dad]; |
|
|
254 first_kid[dad] = v; |
|
|
255 } |
|
|
256 |
|
|
257 /* Depth-first search from dummy root vertex #n */ |
|
|
258 postnum = 0; |
|
|
259 etdfs (n); |
|
|
260 |
|
|
261 SUPERLU_FREE (first_kid); |
|
|
262 SUPERLU_FREE (next_kid); |
|
|
263 return post; |
|
|
264 } |
|
|
265 |
