← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/spanel_dfs.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/spanel_dfs.c
1
 
 
2
 
 
3
 
/*
4
 
 * -- SuperLU routine (version 2.0) --
5
 
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
6
 
 * and Lawrence Berkeley National Lab.
7
 
 * November 15, 1997
8
 
 *
9
 
 */
10
 
/*
11
 
  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
12
 
 
13
 
  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
14
 
  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
15
 
 
16
 
  Permission is hereby granted to use or copy this program for any
17
 
  purpose, provided the above notices are retained on all copies.
18
 
  Permission to modify the code and to distribute modified code is
19
 
  granted, provided the above notices are retained, and a notice that
20
 
  the code was modified is included with the above copyright notice.
21
 
*/
22
 
 
23
 
#include "ssp_defs.h"
24
 
#include "util.h"
25
 
 
26
 
void
27
 
spanel_dfs (
28
 
           const int  m,           /* in - number of rows in the matrix */
29
 
           const int  w,           /* in */
30
 
           const int  jcol,        /* in */
31
 
           SuperMatrix *A,       /* in - original matrix */
32
 
           int        *perm_r,     /* in */
33
 
           int        *nseg,       /* out */
34
 
           float     *dense,      /* out */
35
 
           int        *panel_lsub, /* out */
36
 
           int        *segrep,     /* out */
37
 
           int        *repfnz,     /* out */
38
 
           int        *xprune,     /* out */
39
 
           int        *marker,     /* out */     
40
 
           int        *parent,     /* working array */
41
 
           int        *xplore,     /* working array */
42
 
           GlobalLU_t *Glu         /* modified */
43
 
           )
44
 
{
45
 
/*
46
 
 * Purpose
47
 
 * =======
48
 
 *
49
 
 *   Performs a symbolic factorization on a panel of columns [jcol, jcol+w).
50
 
 *
51
 
 *   A supernode representative is the last column of a supernode.
52
 
 *   The nonzeros in U[*,j] are segments that end at supernodal
53
 
 *   representatives.
54
 
 *
55
 
 *   The routine returns one list of the supernodal representatives
56
 
 *   in topological order of the dfs that generates them. This list is
57
 
 *   a superset of the topological order of each individual column within
58
 
 *   the panel. 
59
 
 *   The location of the first nonzero in each supernodal segment
60
 
 *   (supernodal entry location) is also returned. Each column has a 
61
 
 *   separate list for this purpose.
62
 
 *
63
 
 *   Two marker arrays are used for dfs:
64
 
 *     marker[i] == jj, if i was visited during dfs of current column jj;
65
 
 *     marker1[i] >= jcol, if i was visited by earlier columns in this panel;
66
 
 *
67
 
 *   marker: A-row --> A-row/col (0/1)
68
 
 *   repfnz: SuperA-col --> PA-row
69
 
 *   parent: SuperA-col --> SuperA-col
70
 
 *   xplore: SuperA-col --> index to L-structure
71
 
 *
72
 
 */
73
 
    NCPformat *Astore;
74
 
    float    *a;
75
 
    int       *asub;
76
 
    int       *xa_begin, *xa_end;
77
 
    int       krep, chperm, chmark, chrep, oldrep, kchild, myfnz;
78
 
    int       k, krow, kmark, kperm;
79
 
    int       xdfs, maxdfs, kpar;
80
 
    int       jj;          /* index through each column in the panel */
81
 
    int       *marker1;    /* marker1[jj] >= jcol if vertex jj was visited 
82
 
                              by a previous column within this panel.   */
83
 
    int       *repfnz_col; /* start of each column in the panel */
84
 
    float    *dense_col;  /* start of each column in the panel */
85
 
    int       nextl_col;   /* next available position in panel_lsub[*,jj] */
86
 
    int       *xsup, *supno;
87
 
    int       *lsub, *xlsub;
88
 
 
89
 
    /* Initialize pointers */
90
 
    Astore     = A->Store;
91
 
    a          = Astore->nzval;
92
 
    asub       = Astore->rowind;
93
 
    xa_begin   = Astore->colbeg;
94
 
    xa_end     = Astore->colend;
95
 
    marker1    = marker + m;
96
 
    repfnz_col = repfnz;
97
 
    dense_col  = dense;
98
 
    *nseg      = 0;
99
 
    xsup       = Glu->xsup;
100
 
    supno      = Glu->supno;
101
 
    lsub       = Glu->lsub;
102
 
    xlsub      = Glu->xlsub;
103
 
 
104
 
    /* For each column in the panel */
105
 
    for (jj = jcol; jj < jcol + w; jj++) {
106
 
        nextl_col = (jj - jcol) * m;
107
 
 
108
 
#ifdef CHK_DFS
109
 
        printf("\npanel col %d: ", jj);
110
 
#endif
111
 
 
112
 
        /* For each nonz in A[*,jj] do dfs */
113
 
        for (k = xa_begin[jj]; k < xa_end[jj]; k++) {
114
 
            krow = asub[k];
115
 
            dense_col[krow] = a[k];
116
 
            kmark = marker[krow];       
117
 
            if ( kmark == jj ) 
118
 
                continue;     /* krow visited before, go to the next nonzero */
119
 
 
120
 
            /* For each unmarked nbr krow of jj
121
 
             * krow is in L: place it in structure of L[*,jj]
122
 
             */
123
 
            marker[krow] = jj;
124
 
            kperm = perm_r[krow];
125
 
            
126
 
            if ( kperm == EMPTY ) {
127
 
                panel_lsub[nextl_col++] = krow; /* krow is indexed into A */
128
 
            }
129
 
            /* 
130
 
             * krow is in U: if its supernode-rep krep
131
 
             * has been explored, update repfnz[*]
132
 
             */
133
 
            else {
134
 
                
135
 
                krep = xsup[supno[kperm]+1] - 1;
136
 
                myfnz = repfnz_col[krep];
137
 
                
138
 
#ifdef CHK_DFS
139
 
                printf("krep %d, myfnz %d, perm_r[%d] %d\n", krep, myfnz, krow, kperm);
140
 
#endif
141
 
                if ( myfnz != EMPTY ) { /* Representative visited before */
142
 
                    if ( myfnz > kperm ) repfnz_col[krep] = kperm;
143
 
                    /* continue; */
144
 
                }
145
 
                else {
146
 
                    /* Otherwise, perform dfs starting at krep */
147
 
                    oldrep = EMPTY;
148
 
                    parent[krep] = oldrep;
149
 
                    repfnz_col[krep] = kperm;
150
 
                    xdfs = xlsub[krep];
151
 
                    maxdfs = xprune[krep];
152
 
                    
153
 
#ifdef CHK_DFS 
154
 
                    printf("  xdfs %d, maxdfs %d: ", xdfs, maxdfs);
155
 
                    for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]);
156
 
                    printf("\n");
157
 
#endif
158
 
                    do {
159
 
                        /* 
160
 
                         * For each unmarked kchild of krep 
161
 
                         */
162
 
                        while ( xdfs < maxdfs ) {
163
 
                            
164
 
                            kchild = lsub[xdfs];
165
 
                            xdfs++;
166
 
                            chmark = marker[kchild];
167
 
                            
168
 
                            if ( chmark != jj ) { /* Not reached yet */
169
 
                                marker[kchild] = jj;
170
 
                                chperm = perm_r[kchild];
171
 
                              
172
 
                                /* Case kchild is in L: place it in L[*,j] */
173
 
                                if ( chperm == EMPTY ) {
174
 
                                    panel_lsub[nextl_col++] = kchild;
175
 
                                } 
176
 
                                /* Case kchild is in U: 
177
 
                                 *   chrep = its supernode-rep. If its rep has 
178
 
                                 *   been explored, update its repfnz[*]
179
 
                                 */
180
 
                                else {
181
 
                                    
182
 
                                    chrep = xsup[supno[chperm]+1] - 1;
183
 
                                    myfnz = repfnz_col[chrep];
184
 
#ifdef CHK_DFS
185
 
                                    printf("chrep %d,myfnz %d,perm_r[%d] %d\n",chrep,myfnz,kchild,chperm);
186
 
#endif
187
 
                                    if ( myfnz != EMPTY ) { /* Visited before */
188
 
                                        if ( myfnz > chperm )
189
 
                                            repfnz_col[chrep] = chperm;
190
 
                                    }
191
 
                                    else {
192
 
                                        /* Cont. dfs at snode-rep of kchild */
193
 
                                        xplore[krep] = xdfs;    
194
 
                                        oldrep = krep;
195
 
                                        krep = chrep; /* Go deeper down G(L) */
196
 
                                        parent[krep] = oldrep;
197
 
                                        repfnz_col[krep] = chperm;
198
 
                                        xdfs = xlsub[krep];     
199
 
                                        maxdfs = xprune[krep];
200
 
#ifdef CHK_DFS 
201
 
                                        printf("  xdfs %d, maxdfs %d: ", xdfs, maxdfs);
202
 
                                        for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]); 
203
 
                                        printf("\n");
204
 
#endif
205
 
                                    } /* else */
206
 
                                  
207
 
                                } /* else */
208
 
                              
209
 
                            } /* if... */
210
 
                            
211
 
                        } /* while xdfs < maxdfs */
212
 
                        
213
 
                        /* krow has no more unexplored nbrs:
214
 
                         *    Place snode-rep krep in postorder DFS, if this 
215
 
                         *    segment is seen for the first time. (Note that
216
 
                         *    "repfnz[krep]" may change later.)
217
 
                         *    Backtrack dfs to its parent.
218
 
                         */
219
 
                        if ( marker1[krep] < jcol ) {
220
 
                            segrep[*nseg] = krep;
221
 
                            ++(*nseg);
222
 
                            marker1[krep] = jj;
223
 
                        }
224
 
                        
225
 
                        kpar = parent[krep]; /* Pop stack, mimic recursion */
226
 
                        if ( kpar == EMPTY ) break; /* dfs done */
227
 
                        krep = kpar;
228
 
                        xdfs = xplore[krep];
229
 
                        maxdfs = xprune[krep];
230
 
                        
231
 
#ifdef CHK_DFS 
232
 
                        printf("  pop stack: krep %d,xdfs %d,maxdfs %d: ", krep,xdfs,maxdfs);
233
 
                        for (i = xdfs; i < maxdfs; i++) printf(" %d", lsub[i]);
234
 
                        printf("\n");
235
 
#endif
236
 
                    } while ( kpar != EMPTY ); /* do-while - until empty stack */
237
 
                    
238
 
                } /* else */
239
 
                
240
 
            } /* else */
241
 
            
242
 
        } /* for each nonz in A[*,jj] */
243
 
        
244
 
        repfnz_col += m;    /* Move to next column */
245
 
        dense_col += m;
246
 
        
247
 
    } /* for jj ... */
248
 
    
249
 
}