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