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