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