Line data Source code
1 : /*
2 : Copyright (c) 2003, The Regents of the University of California, through
3 : Lawrence Berkeley National Laboratory (subject to receipt of any required
4 : approvals from U.S. Dept. of Energy)
5 :
6 : All rights reserved.
7 :
8 : The source code is distributed under BSD license, see the file License.txt
9 : at the top-level directory.
10 : */
11 : /*
12 : * -- SuperLU routine (version 3.1) --
13 : * Univ. of California Berkeley, Xerox Palo Alto Research Center,
14 : * and Lawrence Berkeley National Lab.
15 : * August 1, 2008
16 : *
17 : * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
18 : *
19 : * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
20 : * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
21 : *
22 : * Permission is hereby granted to use or copy this program for any
23 : * purpose, provided the above notices are retained on all copies.
24 : * Permission to modify the code and to distribute modified code is
25 : * granted, provided the above notices are retained, and a notice that
26 : * the code was modified is included with the above copyright notice.
27 : */
28 : /*! \file
29 : * \brief Tree layout and computation routines
30 : *
31 : * \ingroup Common
32 : */
33 :
34 : /* Elimination tree computation and layout routines */
35 :
36 : #include <stdio.h>
37 : #include <stdlib.h>
38 : #include "slu_ddefs.h"
39 :
40 : /*
41 : * Implementation of disjoint set union routines.
42 : * Elements are integers in 0..n-1, and the
43 : * names of the sets themselves are of type int.
44 : *
45 : * Calls are:
46 : * initialize_disjoint_sets (n) initial call.
47 : * s = make_set (i) returns a set containing only i.
48 : * s = link (t, u) returns s = t union u, destroying t and u.
49 : * s = find (i) return name of set containing i.
50 : * finalize_disjoint_sets final call.
51 : *
52 : * This implementation uses path compression but not weighted union.
53 : * See Tarjan's book for details.
54 : * John Gilbert, CMI, 1987.
55 : *
56 : * Implemented path-halving by XSL 07/05/95.
57 : */
58 :
59 :
60 : static
61 0 : int *mxCallocInt(int n)
62 : {
63 : register int i;
64 : int *buf;
65 :
66 0 : buf = (int *) SUPERLU_MALLOC( n * sizeof(int) );
67 0 : if ( !buf ) {
68 0 : ABORT("SUPERLU_MALLOC fails for buf in mxCallocInt()");
69 : }
70 0 : for (i = 0; i < n; i++) buf[i] = 0;
71 0 : return (buf);
72 : }
73 :
74 : static
75 : void initialize_disjoint_sets (
76 : int n,
77 : int **pp
78 : )
79 : {
80 0 : (*pp) = mxCallocInt(n);
81 : }
82 :
83 :
84 : static
85 : int make_set (
86 : int i,
87 : int *pp
88 : )
89 : {
90 0 : pp[i] = i;
91 : return i;
92 : }
93 :
94 :
95 : static
96 : int link (
97 : int s,
98 : int t,
99 : int *pp
100 : )
101 : {
102 0 : pp[s] = t;
103 : return t;
104 : }
105 :
106 :
107 : /* PATH HALVING */
108 : static
109 : int find (
110 : int i,
111 : int *pp
112 : )
113 : {
114 : register int p, gp;
115 :
116 0 : p = pp[i];
117 0 : gp = pp[p];
118 0 : while (gp != p) {
119 0 : pp[i] = gp;
120 : i = gp;
121 0 : p = pp[i];
122 0 : gp = pp[p];
123 : }
124 : return (p);
125 : }
126 :
127 : #if 0
128 : /* PATH COMPRESSION */
129 : static
130 : int find (
131 : int i
132 : )
133 : {
134 : if (pp[i] != i)
135 : pp[i] = find (pp[i]);
136 : return pp[i];
137 : }
138 : #endif
139 :
140 : static
141 : void finalize_disjoint_sets (
142 : int *pp
143 : )
144 : {
145 0 : SUPERLU_FREE(pp);
146 : }
147 :
148 :
149 : /*
150 : * Find the elimination tree for A'*A.
151 : * This uses something similar to Liu's algorithm.
152 : * It runs in time O(nz(A)*log n) and does not form A'*A.
153 : *
154 : * Input:
155 : * Sparse matrix A. Numeric values are ignored, so any
156 : * explicit zeros are treated as nonzero.
157 : * Output:
158 : * Integer array of parents representing the elimination
159 : * tree of the symbolic product A'*A. Each vertex is a
160 : * column of A, and nc means a root of the elimination forest.
161 : *
162 : * John R. Gilbert, Xerox, 10 Dec 1990
163 : * Based on code by JRG dated 1987, 1988, and 1990.
164 : */
165 :
166 : /*
167 : * Nonsymmetric elimination tree
168 : */
169 0 : int sp_coletree(
170 : const int_t *acolst, const int_t *acolend, /* column start and end past 1 */
171 : const int_t *arow, /* row indices of A */
172 : int nr, int nc, /* dimension of A */
173 : int *parent) /* parent in elim tree */
174 : {
175 : int *root; /* root of subtree of etree */
176 : int *firstcol; /* first nonzero col in each row*/
177 : int rset, cset;
178 : int row, col;
179 : int rroot;
180 : int p;
181 : int *pp;
182 :
183 0 : root = mxCallocInt (nc);
184 : initialize_disjoint_sets (nc, &pp);
185 :
186 : /* Compute firstcol[row] = first nonzero column in row */
187 :
188 0 : firstcol = mxCallocInt (nr);
189 0 : for (row = 0; row < nr; firstcol[row++] = nc);
190 0 : for (col = 0; col < nc; col++)
191 0 : for (p = acolst[col]; p < acolend[col]; p++) {
192 0 : row = arow[p];
193 0 : firstcol[row] = SUPERLU_MIN(firstcol[row], col);
194 : }
195 :
196 : /* Compute etree by Liu's algorithm for symmetric matrices,
197 : except use (firstcol[r],c) in place of an edge (r,c) of A.
198 : Thus each row clique in A'*A is replaced by a star
199 : centered at its first vertex, which has the same fill. */
200 :
201 0 : for (col = 0; col < nc; col++) {
202 : cset = make_set (col, pp);
203 0 : root[cset] = col;
204 0 : parent[col] = nc; /* Matlab */
205 0 : for (p = acolst[col]; p < acolend[col]; p++) {
206 0 : row = firstcol[arow[p]];
207 0 : if (row >= col) continue;
208 : rset = find (row, pp);
209 0 : rroot = root[rset];
210 0 : if (rroot != col) {
211 0 : parent[rroot] = col;
212 : cset = link (cset, rset, pp);
213 0 : root[cset] = col;
214 : }
215 : }
216 : }
217 :
218 0 : SUPERLU_FREE (root);
219 0 : SUPERLU_FREE (firstcol);
220 : finalize_disjoint_sets (pp);
221 0 : return 0;
222 : }
223 :
224 : /*
225 : * q = TreePostorder (n, p);
226 : *
227 : * Postorder a tree.
228 : * Input:
229 : * p is a vector of parent pointers for a forest whose
230 : * vertices are the integers 0 to n-1; p[root]==n.
231 : * Output:
232 : * q is a vector indexed by 0..n-1 such that q[i] is the
233 : * i-th vertex in a postorder numbering of the tree.
234 : *
235 : * ( 2/7/95 modified by X.Li:
236 : * q is a vector indexed by 0:n-1 such that vertex i is the
237 : * q[i]-th vertex in a postorder numbering of the tree.
238 : * That is, this is the inverse of the previous q. )
239 : *
240 : * In the child structure, lower-numbered children are represented
241 : * first, so that a tree which is already numbered in postorder
242 : * will not have its order changed.
243 : *
244 : * Written by John Gilbert, Xerox, 10 Dec 1990.
245 : * Based on code written by John Gilbert at CMI in 1987.
246 : */
247 :
248 : #if 0 // replaced by a non-recursive version
249 : static
250 : /*
251 : * Depth-first search from vertex v.
252 : */
253 : void etdfs (
254 : int v,
255 : int first_kid[],
256 : int next_kid[],
257 : int post[],
258 : int *postnum
259 : )
260 : {
261 : int w;
262 :
263 : for (w = first_kid[v]; w != -1; w = next_kid[w]) {
264 : etdfs (w, first_kid, next_kid, post, postnum);
265 : }
266 : /* post[postnum++] = v; in Matlab */
267 : post[v] = (*postnum)++; /* Modified by X. Li on 08/10/07 */
268 : }
269 : #endif
270 :
271 : static
272 : /*
273 : * Depth-first search from vertex n. No recursion.
274 : * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
275 : */
276 0 : void nr_etdfs (int n, const int *parent,
277 : const int *first_kid, const int *next_kid,
278 : int *post, int postnum)
279 : {
280 : int current = n, first, next;
281 :
282 0 : while (postnum != n){
283 :
284 : /* no kid for the current node */
285 0 : first = first_kid[current];
286 :
287 : /* no first kid for the current node */
288 0 : if (first == -1){
289 :
290 : /* numbering this node because it has no kid */
291 0 : post[current] = postnum++;
292 :
293 : /* looking for the next kid */
294 0 : next = next_kid[current];
295 :
296 0 : while (next == -1){
297 :
298 : /* no more kids : back to the parent node */
299 0 : current = parent[current];
300 :
301 : /* numbering the parent node */
302 0 : post[current] = postnum++;
303 :
304 : /* get the next kid */
305 0 : next = next_kid[current];
306 : }
307 :
308 : /* stopping criterion */
309 0 : if (postnum==n+1) return;
310 :
311 : /* updating current node */
312 : current = next;
313 : }
314 : /* updating current node */
315 : else {
316 : current = first;
317 : }
318 : }
319 : }
320 :
321 : /*
322 : * Post order a tree
323 : */
324 0 : int *TreePostorder(
325 : int n,
326 : int *parent
327 : )
328 : {
329 : int *first_kid, *next_kid; /* Linked list of children. */
330 : int *post, postnum;
331 : int v, dad;
332 :
333 : /* Allocate storage for working arrays and results */
334 0 : first_kid = mxCallocInt (n+1);
335 0 : next_kid = mxCallocInt (n+1);
336 0 : post = mxCallocInt (n+1);
337 :
338 : /* Set up structure describing children */
339 0 : for (v = 0; v <= n; first_kid[v++] = -1);
340 0 : for (v = n-1; v >= 0; v--) {
341 0 : dad = parent[v];
342 0 : next_kid[v] = first_kid[dad];
343 0 : first_kid[dad] = v;
344 : }
345 :
346 : /* Depth-first search from dummy root vertex #n */
347 : postnum = 0;
348 : #if 0
349 : /* recursion */
350 : etdfs (n, first_kid, next_kid, post, &postnum);
351 : #else
352 : /* no recursion */
353 0 : nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
354 : #endif
355 :
356 0 : SUPERLU_FREE (first_kid);
357 0 : SUPERLU_FREE (next_kid);
358 0 : return post;
359 : }
360 :
361 :
362 : /*
363 : * p = spsymetree (A);
364 : *
365 : * Find the elimination tree for symmetric matrix A.
366 : * This uses Liu's algorithm, and runs in time O(nz*log n).
367 : *
368 : * Input:
369 : * Square sparse matrix A. No check is made for symmetry;
370 : * elements below and on the diagonal are ignored.
371 : * Numeric values are ignored, so any explicit zeros are
372 : * treated as nonzero.
373 : * Output:
374 : * Integer array of parents representing the etree, with n
375 : * meaning a root of the elimination forest.
376 : * Note:
377 : * This routine uses only the upper triangle, while sparse
378 : * Cholesky (as in spchol.c) uses only the lower. Matlab's
379 : * dense Cholesky uses only the upper. This routine could
380 : * be modified to use the lower triangle either by transposing
381 : * the matrix or by traversing it by rows with auxiliary
382 : * pointer and link arrays.
383 : *
384 : * John R. Gilbert, Xerox, 10 Dec 1990
385 : * Based on code by JRG dated 1987, 1988, and 1990.
386 : * Modified by X.S. Li, November 1999.
387 : */
388 :
389 : /*
390 : * Symmetric elimination tree
391 : */
392 0 : int sp_symetree(
393 : const int *acolst, const int *acolend, /* column starts and ends past 1 */
394 : const int *arow, /* row indices of A */
395 : int n, /* dimension of A */
396 : int *parent) /* parent in elim tree */
397 : {
398 : int *root; /* root of subtree of etree */
399 : int rset, cset;
400 : int row, col;
401 : int rroot;
402 : int p;
403 : int *pp;
404 :
405 0 : root = mxCallocInt (n);
406 : initialize_disjoint_sets (n, &pp);
407 :
408 0 : for (col = 0; col < n; col++) {
409 : cset = make_set (col, pp);
410 0 : root[cset] = col;
411 0 : parent[col] = n; /* Matlab */
412 0 : for (p = acolst[col]; p < acolend[col]; p++) {
413 0 : row = arow[p];
414 0 : if (row >= col) continue;
415 : rset = find (row, pp);
416 0 : rroot = root[rset];
417 0 : if (rroot != col) {
418 0 : parent[rroot] = col;
419 : cset = link (cset, rset, pp);
420 0 : root[cset] = col;
421 : }
422 : }
423 : }
424 0 : SUPERLU_FREE (root);
425 : finalize_disjoint_sets (pp);
426 0 : return 0;
427 : } /* SP_SYMETREE */
|
