← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/sp_coletree.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

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