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 ccolumn_dfs.c
13 * \brief Performs a symbolic factorization
14 *
15 * <pre>
16 * -- SuperLU routine (version 3.0) --
17 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
18 * and Lawrence Berkeley National Lab.
19 * October 15, 2003
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 #include "slu_cdefs.h"
35
36 /*! \brief What type of supernodes we want */
37 #define T2_SUPER
38
39
40 /*! \brief
41 *
42 * <pre>
43 * Purpose
44 * =======
45 * CCOLUMN_DFS performs a symbolic factorization on column jcol, and
46 * decide the supernode boundary.
47 *
48 * This routine does not use numeric values, but only use the RHS
49 * row indices to start the dfs.
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. The routine returns a list of such supernodal
54 * representatives in topological order of the dfs that generates them.
55 * The location of the first nonzero in each such supernodal segment
56 * (supernodal entry location) is also returned.
57 *
58 * Local parameters
59 * ================
60 * nseg: no of segments in current U[*,j]
61 * jsuper: jsuper=EMPTY if column j does not belong to the same
62 * supernode as j-1. Otherwise, jsuper=nsuper.
63 *
64 * marker2: A-row --> A-row/col (0/1)
65 * repfnz: SuperA-col --> PA-row
66 * parent: SuperA-col --> SuperA-col
67 * xplore: SuperA-col --> index to L-structure
68 *
69 * Return value
70 * ============
71 * 0 success;
72 * > 0 number of bytes allocated when run out of space.
73 * </pre>
74 */
75 int
ccolumn_dfs(const int m,const int jcol,int * perm_r,int * nseg,int * lsub_col,int * segrep,int * repfnz,int * xprune,int * marker,int * parent,int * xplore,GlobalLU_t * Glu)76 ccolumn_dfs(
77 const int m, /* in - number of rows in the matrix */
78 const int jcol, /* in */
79 int *perm_r, /* in */
80 int *nseg, /* modified - with new segments appended */
81 int *lsub_col, /* in - defines the RHS vector to start the dfs */
82 int *segrep, /* modified - with new segments appended */
83 int *repfnz, /* modified */
84 int *xprune, /* modified */
85 int *marker, /* modified */
86 int *parent, /* working array */
87 int *xplore, /* working array */
88 GlobalLU_t *Glu /* modified */
89 )
90 {
91
92 int jcolp1, jcolm1, jsuper, nsuper, nextl;
93 int k, krep, krow, kmark, kperm;
94 int *marker2; /* Used for small panel LU */
95 int fsupc; /* First column of a snode */
96 int myfnz; /* First nonz column of a U-segment */
97 int chperm, chmark, chrep, kchild;
98 int xdfs, maxdfs, kpar, oldrep;
99 int jptr, jm1ptr;
100 int ito, ifrom, istop; /* Used to compress row subscripts */
101 int mem_error;
102 int *xsup, *supno, *lsub, *xlsub;
103 int nzlmax;
104 int maxsuper;
105
106 xsup = Glu->xsup;
107 supno = Glu->supno;
108 lsub = Glu->lsub;
109 xlsub = Glu->xlsub;
110 nzlmax = Glu->nzlmax;
111
112 maxsuper = sp_ienv(3);
113 jcolp1 = jcol + 1;
114 jcolm1 = jcol - 1;
115 nsuper = supno[jcol];
116 jsuper = nsuper;
117 nextl = xlsub[jcol];
118 marker2 = &marker[2*m];
119
120 /* For each nonzero in A[*,jcol] do dfs */
121 for (k = 0; lsub_col[k] != EMPTY; k++) {
122
123 krow = lsub_col[k];
124 lsub_col[k] = EMPTY;
125 kmark = marker2[krow];
126
127 /* krow was visited before, go to the next nonz */
128 if ( kmark == jcol ) continue;
129
130 /* For each unmarked nbr krow of jcol
131 * krow is in L: place it in structure of L[*,jcol]
132 */
133 marker2[krow] = jcol;
134 kperm = perm_r[krow];
135
136 if ( kperm == EMPTY ) {
137 lsub[nextl++] = krow; /* krow is indexed into A */
138 if ( nextl >= nzlmax ) {
139 if ( mem_error = cLUMemXpand(jcol, nextl, LSUB, &nzlmax, Glu) )
140 return (mem_error);
141 lsub = Glu->lsub;
142 }
143 if ( kmark != jcolm1 ) jsuper = EMPTY;/* Row index subset testing */
144 } else {
145 /* krow is in U: if its supernode-rep krep
146 * has been explored, update repfnz[*]
147 */
148 krep = xsup[supno[kperm]+1] - 1;
149 myfnz = repfnz[krep];
150
151 if ( myfnz != EMPTY ) { /* Visited before */
152 if ( myfnz > kperm ) repfnz[krep] = kperm;
153 /* continue; */
154 }
155 else {
156 /* Otherwise, perform dfs starting at krep */
157 oldrep = EMPTY;
158 parent[krep] = oldrep;
159 repfnz[krep] = kperm;
160 xdfs = xlsub[krep];
161 maxdfs = xprune[krep];
162
163 do {
164 /*
165 * For each unmarked kchild of krep
166 */
167 while ( xdfs < maxdfs ) {
168
169 kchild = lsub[xdfs];
170 xdfs++;
171 chmark = marker2[kchild];
172
173 if ( chmark != jcol ) { /* Not reached yet */
174 marker2[kchild] = jcol;
175 chperm = perm_r[kchild];
176
177 /* Case kchild is in L: place it in L[*,k] */
178 if ( chperm == EMPTY ) {
179 lsub[nextl++] = kchild;
180 if ( nextl >= nzlmax ) {
181 if ( mem_error =
182 cLUMemXpand(jcol,nextl,LSUB,&nzlmax,Glu) )
183 return (mem_error);
184 lsub = Glu->lsub;
185 }
186 if ( chmark != jcolm1 ) jsuper = EMPTY;
187 } else {
188 /* Case kchild is in U:
189 * chrep = its supernode-rep. If its rep has
190 * been explored, update its repfnz[*]
191 */
192 chrep = xsup[supno[chperm]+1] - 1;
193 myfnz = repfnz[chrep];
194 if ( myfnz != EMPTY ) { /* Visited before */
195 if ( myfnz > chperm )
196 repfnz[chrep] = chperm;
197 } else {
198 /* Continue dfs at super-rep of kchild */
199 xplore[krep] = xdfs;
200 oldrep = krep;
201 krep = chrep; /* Go deeper down G(L^t) */
202 parent[krep] = oldrep;
203 repfnz[krep] = chperm;
204 xdfs = xlsub[krep];
205 maxdfs = xprune[krep];
206 } /* else */
207
208 } /* else */
209
210 } /* if */
211
212 } /* while */
213
214 /* krow has no more unexplored nbrs;
215 * place supernode-rep krep in postorder DFS.
216 * backtrack dfs to its parent
217 */
218 segrep[*nseg] = krep;
219 ++(*nseg);
220 kpar = parent[krep]; /* Pop from stack, mimic recursion */
221 if ( kpar == EMPTY ) break; /* dfs done */
222 krep = kpar;
223 xdfs = xplore[krep];
224 maxdfs = xprune[krep];
225
226 } while ( kpar != EMPTY ); /* Until empty stack */
227
228 } /* else */
229
230 } /* else */
231
232 } /* for each nonzero ... */
233
234 /* Check to see if j belongs in the same supernode as j-1 */
235 if ( jcol == 0 ) { /* Do nothing for column 0 */
236 nsuper = supno[0] = 0;
237 } else {
238 fsupc = xsup[nsuper];
239 jptr = xlsub[jcol]; /* Not compressed yet */
240 jm1ptr = xlsub[jcolm1];
241
242 #ifdef T2_SUPER
243 if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = EMPTY;
244 #endif
245 /* Make sure the number of columns in a supernode doesn't
246 exceed threshold. */
247 if ( jcol - fsupc >= maxsuper ) jsuper = EMPTY;
248
249 /* If jcol starts a new supernode, reclaim storage space in
250 * lsub from the previous supernode. Note we only store
251 * the subscript set of the first and last columns of
252 * a supernode. (first for num values, last for pruning)
253 */
254 if ( jsuper == EMPTY ) { /* starts a new supernode */
255 if ( (fsupc < jcolm1-1) ) { /* >= 3 columns in nsuper */
256 #ifdef CHK_COMPRESS
257 printf(" Compress lsub[] at super %d-%d\n", fsupc, jcolm1);
258 #endif
259 ito = xlsub[fsupc+1];
260 xlsub[jcolm1] = ito;
261 istop = ito + jptr - jm1ptr;
262 xprune[jcolm1] = istop; /* Initialize xprune[jcol-1] */
263 xlsub[jcol] = istop;
264 for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
265 lsub[ito] = lsub[ifrom];
266 nextl = ito; /* = istop + length(jcol) */
267 }
268 nsuper++;
269 supno[jcol] = nsuper;
270 } /* if a new supernode */
271
272 } /* else: jcol > 0 */
273
274 /* Tidy up the pointers before exit */
275 xsup[nsuper+1] = jcolp1;
276 supno[jcolp1] = nsuper;
277 xprune[jcol] = nextl; /* Initialize upper bound for pruning */
278 xlsub[jcolp1] = nextl;
279
280 return 0;
281 }
282