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 dcolumn_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_ddefs.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 : * DCOLUMN_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=SLU_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
76 0 : dcolumn_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_t *xprune, /* modified */
85 : int *marker, /* modified */
86 : int *parent, /* working array */
87 : int_t *xplore, /* working array */
88 : GlobalLU_t *Glu /* modified */
89 : )
90 : {
91 :
92 : int jcolp1, jcolm1, jsuper, nsuper;
93 : int 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_t xdfs, maxdfs, nextl, k;
99 : int kpar, oldrep;
100 : int_t jptr, jm1ptr;
101 : int_t ito, ifrom, istop; /* Used to compress row subscripts */
102 : int *xsup, *supno;
103 : int_t *lsub, *xlsub;
104 : int_t nzlmax, mem_error;
105 : int maxsuper;
106 :
107 0 : xsup = Glu->xsup;
108 0 : supno = Glu->supno;
109 0 : lsub = Glu->lsub;
110 0 : xlsub = Glu->xlsub;
111 0 : nzlmax = Glu->nzlmax;
112 :
113 0 : maxsuper = sp_ienv(3);
114 0 : jcolp1 = jcol + 1;
115 0 : jcolm1 = jcol - 1;
116 0 : nsuper = supno[jcol];
117 : jsuper = nsuper;
118 0 : nextl = xlsub[jcol];
119 0 : marker2 = &marker[2*m];
120 :
121 : /* For each nonzero in A[*,jcol] do dfs */
122 0 : for (k = 0; lsub_col[k] != SLU_EMPTY; k++) {
123 :
124 : krow = lsub_col[k];
125 0 : lsub_col[k] = SLU_EMPTY;
126 0 : kmark = marker2[krow];
127 :
128 : /* krow was visited before, go to the next nonz */
129 0 : if ( kmark == jcol ) continue;
130 :
131 : /* For each unmarked nbr krow of jcol
132 : * krow is in L: place it in structure of L[*,jcol]
133 : */
134 0 : marker2[krow] = jcol;
135 0 : kperm = perm_r[krow];
136 :
137 0 : if ( kperm == SLU_EMPTY ) {
138 0 : lsub[nextl++] = krow; /* krow is indexed into A */
139 0 : if ( nextl >= nzlmax ) {
140 0 : mem_error = dLUMemXpand(jcol, nextl, LSUB, &nzlmax, Glu);
141 0 : if ( mem_error ) return (mem_error);
142 0 : lsub = Glu->lsub;
143 : }
144 0 : if ( kmark != jcolm1 ) jsuper = SLU_EMPTY;/* Row index subset testing */
145 : } else {
146 : /* krow is in U: if its supernode-rep krep
147 : * has been explored, update repfnz[*]
148 : */
149 0 : krep = xsup[supno[kperm]+1] - 1;
150 0 : myfnz = repfnz[krep];
151 :
152 0 : if ( myfnz != SLU_EMPTY ) { /* Visited before */
153 0 : if ( myfnz > kperm ) repfnz[krep] = kperm;
154 : /* continue; */
155 : }
156 : else {
157 : /* Otherwise, perform dfs starting at krep */
158 : oldrep = SLU_EMPTY;
159 0 : parent[krep] = oldrep;
160 0 : repfnz[krep] = kperm;
161 0 : xdfs = xlsub[krep];
162 0 : maxdfs = xprune[krep];
163 :
164 : do {
165 : /*
166 : * For each unmarked kchild of krep
167 : */
168 0 : while ( xdfs < maxdfs ) {
169 :
170 0 : kchild = lsub[xdfs];
171 0 : xdfs++;
172 0 : chmark = marker2[kchild];
173 :
174 0 : if ( chmark != jcol ) { /* Not reached yet */
175 0 : marker2[kchild] = jcol;
176 0 : chperm = perm_r[kchild];
177 :
178 : /* Case kchild is in L: place it in L[*,k] */
179 0 : if ( chperm == SLU_EMPTY ) {
180 0 : lsub[nextl++] = kchild;
181 0 : if ( nextl >= nzlmax ) {
182 : mem_error =
183 0 : dLUMemXpand(jcol,nextl,LSUB,&nzlmax,Glu);
184 0 : if ( mem_error ) return (mem_error);
185 0 : lsub = Glu->lsub;
186 : }
187 0 : if ( chmark != jcolm1 ) jsuper = SLU_EMPTY;
188 : } else {
189 : /* Case kchild is in U:
190 : * chrep = its supernode-rep. If its rep has
191 : * been explored, update its repfnz[*]
192 : */
193 0 : chrep = xsup[supno[chperm]+1] - 1;
194 0 : myfnz = repfnz[chrep];
195 0 : if ( myfnz != SLU_EMPTY ) { /* Visited before */
196 0 : if ( myfnz > chperm )
197 0 : repfnz[chrep] = chperm;
198 : } else {
199 : /* Continue dfs at super-rep of kchild */
200 0 : xplore[krep] = xdfs;
201 : oldrep = krep;
202 : krep = chrep; /* Go deeper down G(L^t) */
203 0 : parent[krep] = oldrep;
204 0 : repfnz[krep] = chperm;
205 0 : xdfs = xlsub[krep];
206 0 : maxdfs = xprune[krep];
207 : } /* else */
208 :
209 : } /* else */
210 :
211 : } /* if */
212 :
213 : } /* while */
214 :
215 : /* krow has no more unexplored nbrs;
216 : * place supernode-rep krep in postorder DFS.
217 : * backtrack dfs to its parent
218 : */
219 0 : segrep[*nseg] = krep;
220 0 : ++(*nseg);
221 0 : kpar = parent[krep]; /* Pop from stack, mimic recursion */
222 0 : if ( kpar == SLU_EMPTY ) break; /* dfs done */
223 : krep = kpar;
224 0 : xdfs = xplore[krep];
225 0 : maxdfs = xprune[krep];
226 :
227 : } while ( kpar != SLU_EMPTY ); /* Until empty stack */
228 :
229 : } /* else */
230 :
231 : } /* else */
232 :
233 : } /* for each nonzero ... */
234 :
235 : /* Check to see if j belongs in the same supernode as j-1 */
236 0 : if ( jcol == 0 ) { /* Do nothing for column 0 */
237 0 : nsuper = supno[0] = 0;
238 : } else {
239 0 : fsupc = xsup[nsuper];
240 0 : jptr = xlsub[jcol]; /* Not compressed yet */
241 0 : jm1ptr = xlsub[jcolm1];
242 :
243 : #ifdef T2_SUPER
244 0 : if ( (nextl-jptr != jptr-jm1ptr-1) ) jsuper = SLU_EMPTY;
245 : #endif
246 : /* Make sure the number of columns in a supernode doesn't
247 : exceed threshold. */
248 0 : if ( jcol - fsupc >= maxsuper ) jsuper = SLU_EMPTY;
249 :
250 : /* If jcol starts a new supernode, reclaim storage space in
251 : * lsub from the previous supernode. Note we only store
252 : * the subscript set of the first and last columns of
253 : * a supernode. (first for num values, last for pruning)
254 : */
255 0 : if ( jsuper == SLU_EMPTY ) { /* starts a new supernode */
256 0 : if ( (fsupc < jcolm1-1) ) { /* >= 3 columns in nsuper */
257 : #ifdef CHK_COMPRESS
258 : printf(" Compress lsub[] at super %d-%d\n", fsupc, jcolm1);
259 : #endif
260 0 : ito = xlsub[fsupc+1];
261 0 : xlsub[jcolm1] = ito;
262 0 : istop = ito + jptr - jm1ptr;
263 0 : xprune[jcolm1] = istop; /* Initialize xprune[jcol-1] */
264 0 : xlsub[jcol] = istop;
265 0 : for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
266 0 : lsub[ito] = lsub[ifrom];
267 : nextl = ito; /* = istop + length(jcol) */
268 : }
269 0 : nsuper++;
270 0 : supno[jcol] = nsuper;
271 : } /* if a new supernode */
272 :
273 : } /* else: jcol > 0 */
274 :
275 : /* Tidy up the pointers before exit */
276 0 : xsup[nsuper+1] = jcolp1;
277 0 : supno[jcolp1] = nsuper;
278 0 : xprune[jcol] = nextl; /* Initialize upper bound for pruning */
279 0 : xlsub[jcolp1] = nextl;
280 :
281 0 : return 0;
282 : }
|
