comparison main/sparse/SuperLU/SRC/zpanel_dfs.c @ 0:6b33357c7561 octave-forge

Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
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 "zsp_defs.h"
24 #include "util.h"
25
26 void
27 zpanel_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 doublecomplex *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 doublecomplex *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 doublecomplex *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 }