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