~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
ImportĀ upstreamĀ versionĀ 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
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 "csp_defs.h"
 
24
#include "util.h"
 
25
 
 
26
void
 
27
cpanel_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
           complex     *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
    complex    *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
    complex    *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
}