|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === UMF_analyze ========================================================== */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* -------------------------------------------------------------------------- */ |
|
|
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */ |
|
|
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */ |
|
|
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */ |
|
|
9 /* -------------------------------------------------------------------------- */ |
|
|
10 |
|
|
11 /* |
|
|
12 Symbolic LL' factorization of A'*A, to get upper bounds on the size of |
|
|
13 L and U for LU = PAQ, and to determine the frontal matrices and |
|
|
14 (supernodal) column elimination tree. No fill-reducing column pre-ordering |
|
|
15 is used. |
|
|
16 |
|
|
17 Returns TRUE if successful, FALSE if out of memory. UMF_analyze can only |
|
|
18 run out of memory if anzmax (which is Ap [n_row]) is too small. |
|
|
19 |
|
|
20 Uses workspace of size O(nonzeros in A). On input, the matrix A is |
|
|
21 stored in row-form at the tail end of Ai. It is destroyed on output. |
|
|
22 The rows of A must be sorted by increasing first column index. |
|
|
23 The matrix is assumed to be valid. |
|
|
24 |
|
|
25 Empty rows and columns have already been removed. |
|
|
26 |
|
|
27 */ |
|
|
28 |
|
|
29 #include "umf_internal.h" |
|
|
30 #include "umf_apply_order.h" |
|
|
31 #include "umf_fsize.h" |
|
|
32 |
|
|
33 /* ========================================================================== */ |
|
|
34 |
|
|
35 GLOBAL Int UMF_analyze |
|
|
36 ( |
|
|
37 Int n_row, /* A is n_row-by-n_col */ |
|
|
38 Int n_col, |
|
|
39 Int Ai [ ], /* Ai [Ap [0]..Ap[n_row]-1]: column indices */ |
|
|
40 /* destroyed on output. Note that this is NOT the */ |
|
|
41 /* user's Ai that was passed to UMFPACK_*symbolic */ |
|
|
42 /* size of Ai, Ap [n_row] = anzmax >= anz + n_col */ |
|
|
43 /* Ap [0] must be => n_col. The space to the */ |
|
|
44 /* front of Ai is used as workspace. */ |
|
|
45 |
|
|
46 Int Ap [ ], /* of size MAX (n_row, n_col) + 1 */ |
|
|
47 /* Ap [0..n_row]: row pointers */ |
|
|
48 /* Row i is in Ai [Ap [i] ... Ap [i+1]-1] */ |
|
|
49 |
|
|
50 /* rows must have smallest col index first, or be */ |
|
|
51 /* in sorted form. Used as workspace of size n_col */ |
|
|
52 /* and destroyed. */ |
|
|
53 |
|
|
54 /* Note that this is NOT the */ |
|
|
55 /* user's Ap that was passed to UMFPACK_*symbolic */ |
|
|
56 |
|
|
57 Int Up [ ], /* workspace of size n_col, and output column perm. |
|
|
58 * for column etree postorder. */ |
|
|
59 |
|
|
60 Int fixQ, |
|
|
61 |
|
|
62 /* temporary workspaces: */ |
|
|
63 Int W [ ], /* W [0..n_col-1] */ |
|
|
64 Int Link [ ], /* Link [0..n_col-1] */ |
|
|
65 |
|
|
66 /* output: information about each frontal matrix: */ |
|
|
67 Int Front_ncols [ ], /* size n_col */ |
|
|
68 Int Front_nrows [ ], /* of size n_col */ |
|
|
69 Int Front_npivcol [ ], /* of size n_col */ |
|
|
70 Int Front_parent [ ], /* of size n_col */ |
|
|
71 Int *nfr_out, |
|
|
72 |
|
|
73 Int *p_ncompactions /* number of compactions in UMF_analyze */ |
|
|
74 ) |
|
|
75 { |
|
|
76 /* ====================================================================== */ |
|
|
77 /* ==== local variables ================================================= */ |
|
|
78 /* ====================================================================== */ |
|
|
79 |
|
|
80 Int j, j3, col, k, row, parent, j2, pdest, p, p2, thickness, npivots, nfr, |
|
|
81 i, *Winv, kk, npiv, jnext, krow, knext, pfirst, jlast, ncompactions, |
|
|
82 *Front_stack, *Front_order, *Front_child, *Front_sibling, |
|
|
83 Wflag, npivcol, fallrows, fallcols, fpiv, frows, fcols, *Front_size ; |
|
|
84 |
|
|
85 nfr = 0 ; |
|
|
86 DEBUG0 (("UMF_analyze: anzmax "ID" anrow "ID" ancol "ID"\n", |
|
|
87 Ap [n_row], n_row, n_col)) ; |
|
|
88 |
|
|
89 /* ====================================================================== */ |
|
|
90 /* ==== initializations ================================================= */ |
|
|
91 /* ====================================================================== */ |
|
|
92 |
|
|
93 #pragma ivdep |
|
|
94 for (j = 0 ; j < n_col ; j++) |
|
|
95 { |
|
|
96 Link [j] = EMPTY ; |
|
|
97 W [j] = EMPTY ; |
|
|
98 Up [j] = EMPTY ; |
|
|
99 |
|
|
100 /* Frontal matrix data structure: */ |
|
|
101 Front_npivcol [j] = 0 ; /* number of pivot columns */ |
|
|
102 Front_nrows [j] = 0 ; /* number of rows, incl. pivot rows */ |
|
|
103 Front_ncols [j] = 0 ; /* number of cols, incl. pivot cols */ |
|
|
104 Front_parent [j] = EMPTY ; /* parent front */ |
|
|
105 /* Note that only non-pivotal columns are stored in a front (a "row" */ |
|
|
106 /* of U) during elimination. */ |
|
|
107 } |
|
|
108 |
|
|
109 /* the rows must be sorted by increasing min col */ |
|
|
110 krow = 0 ; |
|
|
111 pfirst = Ap [0] ; |
|
|
112 jlast = EMPTY ; |
|
|
113 jnext = EMPTY ; |
|
|
114 Wflag = 0 ; |
|
|
115 |
|
|
116 /* this test requires the size of Ai to be >= n_col + nz */ |
|
|
117 ASSERT (pfirst >= n_col) ; /* Ai must be large enough */ |
|
|
118 |
|
|
119 /* pdest points to the first free space in Ai */ |
|
|
120 pdest = 0 ; |
|
|
121 ncompactions = 0 ; |
|
|
122 |
|
|
123 /* ====================================================================== */ |
|
|
124 /* === compute symbolic LL' factorization (unsorted) ==================== */ |
|
|
125 /* ====================================================================== */ |
|
|
126 |
|
|
127 for (j = 0 ; j < n_col ; j = jnext) |
|
|
128 { |
|
|
129 DEBUG1 (("\n\n============Front "ID" starting. nfr = "ID"\n", j, nfr)) ; |
|
|
130 |
|
|
131 /* ================================================================== */ |
|
|
132 /* === garbage collection =========================================== */ |
|
|
133 /* ================================================================== */ |
|
|
134 |
|
|
135 if (pdest + (n_col-j) > pfirst) |
|
|
136 { |
|
|
137 /* we might run out ... compact the rows of U */ |
|
|
138 |
|
|
139 #ifndef NDEBUG |
|
|
140 DEBUG0 (("UMF_analyze COMPACTION, j="ID" pfirst="ID"\n", |
|
|
141 j, pfirst)) ; |
|
|
142 for (row = 0 ; row < j ; row++) |
|
|
143 { |
|
|
144 if (Up [row] != EMPTY) |
|
|
145 { |
|
|
146 /* this is a live row of U */ |
|
|
147 DEBUG1 (("Live row: "ID" cols: ", row)) ; |
|
|
148 p = Up [row] ; |
|
|
149 ASSERT (Front_ncols [row] > Front_npivcol [row]) ; |
|
|
150 p2 = p + (Front_ncols [row] - Front_npivcol [row]) ; |
|
|
151 for ( ; p < p2 ; p++) |
|
|
152 { |
|
|
153 DEBUG1 ((ID, Ai [p])) ; |
|
|
154 ASSERT (p < pfirst) ; |
|
|
155 ASSERT (Ai [p] > row && Ai [p] < n_col) ; |
|
|
156 } |
|
|
157 DEBUG1 (("\n")) ; |
|
|
158 } |
|
|
159 } |
|
|
160 DEBUG1 (("\nStarting to compact:\n")) ; |
|
|
161 #endif |
|
|
162 |
|
|
163 pdest = 0 ; |
|
|
164 ncompactions++ ; |
|
|
165 for (row = 0 ; row < j ; row++) |
|
|
166 { |
|
|
167 if (Up [row] != EMPTY) |
|
|
168 { |
|
|
169 /* this is a live row of U */ |
|
|
170 DEBUG1 (("Live row: "ID" cols: ", row)) ; |
|
|
171 ASSERT (row < n_col) ; |
|
|
172 p = Up [row] ; |
|
|
173 ASSERT (Front_ncols [row] > Front_npivcol [row]) ; |
|
|
174 p2 = p + (Front_ncols [row] - Front_npivcol [row]) ; |
|
|
175 Up [row] = pdest ; |
|
|
176 for ( ; p < p2 ; p++) |
|
|
177 { |
|
|
178 DEBUG1 ((ID, Ai [p])) ; |
|
|
179 ASSERT (p < pfirst) ; |
|
|
180 ASSERT (Ai [p] > row && Ai [p] < n_col) ; |
|
|
181 Ai [pdest++] = Ai [p] ; |
|
|
182 ASSERT (pdest <= pfirst) ; |
|
|
183 } |
|
|
184 DEBUG1 (("\n")) ; |
|
|
185 } |
|
|
186 } |
|
|
187 |
|
|
188 #ifndef NDEBUG |
|
|
189 DEBUG1 (("\nAFTER COMPACTION, j="ID" pfirst="ID"\n", j, pfirst)) ; |
|
|
190 for (row = 0 ; row < j ; row++) |
|
|
191 { |
|
|
192 if (Up [row] != EMPTY) |
|
|
193 { |
|
|
194 /* this is a live row of U */ |
|
|
195 DEBUG1 (("Live row: "ID" cols: ", row)) ; |
|
|
196 p = Up [row] ; |
|
|
197 ASSERT (Front_ncols [row] > Front_npivcol [row]) ; |
|
|
198 p2 = p + (Front_ncols [row] - Front_npivcol [row]) ; |
|
|
199 for ( ; p < p2 ; p++) |
|
|
200 { |
|
|
201 DEBUG1 ((ID, Ai [p])) ; |
|
|
202 ASSERT (p < pfirst) ; |
|
|
203 ASSERT (Ai [p] > row && Ai [p] < n_col) ; |
|
|
204 } |
|
|
205 DEBUG1 (("\n")) ; |
|
|
206 } |
|
|
207 } |
|
|
208 #endif |
|
|
209 |
|
|
210 } |
|
|
211 |
|
|
212 if (pdest + (n_col-j) > pfirst) |
|
|
213 { |
|
|
214 /* :: out of memory in umf_analyze :: */ |
|
|
215 /* it can't happen, if pfirst >= n_col */ |
|
|
216 return (FALSE) ; /* internal error! */ |
|
|
217 } |
|
|
218 |
|
|
219 /* ------------------------------------------------------------------ */ |
|
|
220 /* is the last front a child of this one? */ |
|
|
221 /* ------------------------------------------------------------------ */ |
|
|
222 |
|
|
223 if (jlast != EMPTY && Link [j] == jlast) |
|
|
224 { |
|
|
225 /* yes - create row j by appending to jlast */ |
|
|
226 DEBUG1 (("GOT:last front is child of this one: j "ID" jlast "ID"\n", |
|
|
227 j, jlast)) ; |
|
|
228 ASSERT (jlast >= 0 && jlast < j) ; |
|
|
229 |
|
|
230 Up [j] = Up [jlast] ; |
|
|
231 Up [jlast] = EMPTY ; |
|
|
232 |
|
|
233 /* find the parent, delete column j, and update W */ |
|
|
234 parent = n_col ; |
|
|
235 for (p = Up [j] ; p < pdest ; ) |
|
|
236 { |
|
|
237 j3 = Ai [p] ; |
|
|
238 DEBUG1 (("Initial row of U: col "ID" ", j3)) ; |
|
|
239 ASSERT (j3 >= 0 && j3 < n_col) ; |
|
|
240 DEBUG1 (("W: "ID" \n", W [j3])) ; |
|
|
241 ASSERT (W [j3] == Wflag) ; |
|
|
242 if (j == j3) |
|
|
243 { |
|
|
244 DEBUG1 (("Found column j at p = "ID"\n", p)) ; |
|
|
245 Ai [p] = Ai [--pdest] ; |
|
|
246 } |
|
|
247 else |
|
|
248 { |
|
|
249 if (j3 < parent) |
|
|
250 { |
|
|
251 parent = j3 ; |
|
|
252 } |
|
|
253 p++ ; |
|
|
254 } |
|
|
255 } |
|
|
256 |
|
|
257 /* delete jlast from the link list of j */ |
|
|
258 Link [j] = Link [jlast] ; |
|
|
259 |
|
|
260 ASSERT (Front_nrows [jlast] > Front_npivcol [jlast]) ; |
|
|
261 thickness = (Front_nrows [jlast] - Front_npivcol [jlast]) ; |
|
|
262 DEBUG1 (("initial thickness: "ID"\n", thickness)) ; |
|
|
263 |
|
|
264 } |
|
|
265 else |
|
|
266 { |
|
|
267 Up [j] = pdest ; |
|
|
268 parent = n_col ; |
|
|
269 /* thickness: number of (nonpivotal) rows in frontal matrix j */ |
|
|
270 thickness = 0 ; |
|
|
271 Wflag = j ; |
|
|
272 } |
|
|
273 |
|
|
274 /* ================================================================== */ |
|
|
275 /* === compute row j of A*A' ======================================== */ |
|
|
276 /* ================================================================== */ |
|
|
277 |
|
|
278 /* ------------------------------------------------------------------ */ |
|
|
279 /* flag the diagonal entry in row U, but do not add to pattern */ |
|
|
280 /* ------------------------------------------------------------------ */ |
|
|
281 |
|
|
282 ASSERT (pdest <= pfirst) ; |
|
|
283 W [j] = Wflag ; |
|
|
284 |
|
|
285 DEBUG1 (("\nComputing row "ID" of A'*A\n", j)) ; |
|
|
286 DEBUG2 ((" col: "ID" (diagonal)\n", j)) ; |
|
|
287 |
|
|
288 /* ------------------------------------------------------------------ */ |
|
|
289 /* find the rows the contribute to this column j */ |
|
|
290 /* ------------------------------------------------------------------ */ |
|
|
291 |
|
|
292 jnext = n_col ; |
|
|
293 for (knext = krow ; knext < n_row ; knext++) |
|
|
294 { |
|
|
295 ASSERT (Ap [knext] < Ap [knext+1]) ; |
|
|
296 ASSERT (Ap [knext] >= pfirst && Ap [knext] <= Ap [n_row]) ; |
|
|
297 jnext = Ai [Ap [knext]] ; |
|
|
298 ASSERT (jnext >= j) ; |
|
|
299 if (jnext != j) |
|
|
300 { |
|
|
301 break ; |
|
|
302 } |
|
|
303 } |
|
|
304 |
|
|
305 /* rows krow ... knext-1 all have first column index of j */ |
|
|
306 /* (or are empty) */ |
|
|
307 |
|
|
308 /* row knext has first column index of jnext */ |
|
|
309 /* if knext = n_row, then jnext is n_col */ |
|
|
310 if (knext == n_row) |
|
|
311 { |
|
|
312 jnext = n_col ; |
|
|
313 } |
|
|
314 |
|
|
315 ASSERT (jnext > j) ; |
|
|
316 ASSERT (jnext <= n_col) ; |
|
|
317 |
|
|
318 /* ------------------------------------------------------------------ */ |
|
|
319 /* for each nonzero A (k,j) in column j of A do: */ |
|
|
320 /* ------------------------------------------------------------------ */ |
|
|
321 |
|
|
322 for (k = krow ; k < knext ; k++) |
|
|
323 { |
|
|
324 p = Ap [k] ; |
|
|
325 p2 = Ap [k+1] ; |
|
|
326 ASSERT (p < p2) ; |
|
|
327 |
|
|
328 /* merge row k of A into W */ |
|
|
329 DEBUG2 ((" ---- A row "ID" ", k)) ; |
|
|
330 ASSERT (k >= 0 && k < n_row) ; |
|
|
331 ASSERT (Ai [p] == j) ; |
|
|
332 DEBUG2 ((" p "ID" p2 "ID"\n cols:", p, p2)) ; |
|
|
333 ASSERT (p >= pfirst && p < Ap [n_row]) ; |
|
|
334 ASSERT (p2 > pfirst && p2 <= Ap [n_row]) ; |
|
|
335 for ( ; p < p2 ; p++) |
|
|
336 { |
|
|
337 /* add to pattern if seen for the first time */ |
|
|
338 col = Ai [p] ; |
|
|
339 ASSERT (col >= j && col < n_col) ; |
|
|
340 DEBUG3 ((" "ID, col)) ; |
|
|
341 if (W [col] != Wflag) |
|
|
342 { |
|
|
343 Ai [pdest++] = col ; |
|
|
344 ASSERT (pdest <= pfirst) ; |
|
|
345 /* flag this column has having been seen for row j */ |
|
|
346 W [col] = Wflag ; |
|
|
347 if (col < parent) |
|
|
348 { |
|
|
349 parent = col ; |
|
|
350 } |
|
|
351 } |
|
|
352 } |
|
|
353 DEBUG2 (("\n")) ; |
|
|
354 thickness++ ; |
|
|
355 } |
|
|
356 |
|
|
357 #ifndef NDEBUG |
|
|
358 DEBUG3 (("\nRow "ID" of A'A:\n", j)) ; |
|
|
359 for (p = Up [j] ; p < pdest ; p++) |
|
|
360 { |
|
|
361 DEBUG3 ((" "ID, Ai [p])) ; |
|
|
362 } |
|
|
363 DEBUG3 (("\n")) ; |
|
|
364 #endif |
|
|
365 |
|
|
366 /* ------------------------------------------------------------------ */ |
|
|
367 /* delete rows up to but not including knext */ |
|
|
368 /* ------------------------------------------------------------------ */ |
|
|
369 |
|
|
370 krow = knext ; |
|
|
371 pfirst = Ap [knext] ; |
|
|
372 |
|
|
373 /* we can now use Ai [0..pfirst-1] as workspace for rows of U */ |
|
|
374 |
|
|
375 /* ================================================================== */ |
|
|
376 /* === compute jth row of U ========================================= */ |
|
|
377 /* ================================================================== */ |
|
|
378 |
|
|
379 /* for each nonzero U (k,j) in column j of U (1:j-1,:) do */ |
|
|
380 for (k = Link [j] ; k != EMPTY ; k = Link [k]) |
|
|
381 { |
|
|
382 /* merge row k of U into W */ |
|
|
383 DEBUG2 ((" ---- U row "ID, k)) ; |
|
|
384 ASSERT (k >= 0 && k < n_col) ; |
|
|
385 ASSERT (Up [k] != EMPTY) ; |
|
|
386 p = Up [k] ; |
|
|
387 ASSERT (Front_ncols [k] > Front_npivcol [k]) ; |
|
|
388 p2 = p + (Front_ncols [k] - Front_npivcol [k]) ; |
|
|
389 DEBUG2 ((" p "ID" p2 "ID"\n cols:", p, p2)) ; |
|
|
390 ASSERT (p <= pfirst) ; |
|
|
391 ASSERT (p2 <= pfirst) ; |
|
|
392 for ( ; p < p2 ; p++) |
|
|
393 { |
|
|
394 /* add to pattern if seen for the first time */ |
|
|
395 col = Ai [p] ; |
|
|
396 ASSERT (col >= j && col < n_col) ; |
|
|
397 DEBUG3 ((" "ID, col)) ; |
|
|
398 if (W [col] != Wflag) |
|
|
399 { |
|
|
400 Ai [pdest++] = col ; |
|
|
401 ASSERT (pdest <= pfirst) ; |
|
|
402 /* flag this col has having been seen for row j */ |
|
|
403 W [col] = Wflag ; |
|
|
404 if (col < parent) |
|
|
405 { |
|
|
406 parent = col ; |
|
|
407 } |
|
|
408 } |
|
|
409 } |
|
|
410 DEBUG2 (("\n")) ; |
|
|
411 |
|
|
412 /* mark the row k as deleted */ |
|
|
413 Up [k] = EMPTY ; |
|
|
414 |
|
|
415 ASSERT (Front_nrows [k] > Front_npivcol [k]) ; |
|
|
416 thickness += (Front_nrows [k] - Front_npivcol [k]) ; |
|
|
417 ASSERT (Front_parent [k] == j) ; |
|
|
418 } |
|
|
419 |
|
|
420 #ifndef NDEBUG |
|
|
421 DEBUG3 (("\nRow "ID" of U prior to supercolumn detection:\n", j)); |
|
|
422 for (p = Up [j] ; p < pdest ; p++) |
|
|
423 { |
|
|
424 DEBUG3 ((" "ID, Ai [p])) ; |
|
|
425 } |
|
|
426 DEBUG3 (("\n")) ; |
|
|
427 DEBUG1 (("thickness, prior to supercol detect: "ID"\n", thickness)) ; |
|
|
428 #endif |
|
|
429 |
|
|
430 /* ================================================================== */ |
|
|
431 /* === quicky mass elimination ====================================== */ |
|
|
432 /* ================================================================== */ |
|
|
433 |
|
|
434 /* this code detects some supernodes, but it might miss */ |
|
|
435 /* some because the elimination tree (created on the fly) */ |
|
|
436 /* is not yet post-ordered, and because the pattern of A'*A */ |
|
|
437 /* is also computed on the fly. */ |
|
|
438 |
|
|
439 /* j2 is incremented because the pivot columns are not stored */ |
|
|
440 |
|
|
441 for (j2 = j+1 ; j2 < jnext ; j2++) |
|
|
442 { |
|
|
443 ASSERT (j2 >= 0 && j2 < n_col) ; |
|
|
444 if (W [j2] != Wflag || Link [j2] != EMPTY) |
|
|
445 { |
|
|
446 break ; |
|
|
447 } |
|
|
448 } |
|
|
449 |
|
|
450 /* the loop above terminated with j2 at the first non-supernode */ |
|
|
451 DEBUG1 (("jnext = "ID"\n", jnext)) ; |
|
|
452 ASSERT (j2 <= jnext) ; |
|
|
453 jnext = j2 ; |
|
|
454 j2-- ; |
|
|
455 DEBUG1 (("j2 = "ID"\n", j2)) ; |
|
|
456 ASSERT (j2 < n_col) ; |
|
|
457 |
|
|
458 npivots = j2-j+1 ; |
|
|
459 DEBUG1 (("Number of pivot columns: "ID"\n", npivots)) ; |
|
|
460 |
|
|
461 /* rows j:j2 have the same nonzero pattern, except for columns j:j2-1 */ |
|
|
462 |
|
|
463 if (j2 > j) |
|
|
464 { |
|
|
465 /* supernode detected, prune the pattern of new row j */ |
|
|
466 ASSERT (parent == j+1) ; |
|
|
467 ASSERT (j2 < n_col) ; |
|
|
468 DEBUG1 (("Supernode detected, j "ID" to j2 "ID"\n", j, j2)) ; |
|
|
469 |
|
|
470 parent = n_col ; |
|
|
471 p2 = pdest ; |
|
|
472 pdest = Up [j] ; |
|
|
473 for (p = Up [j] ; p < p2 ; p++) |
|
|
474 { |
|
|
475 col = Ai [p] ; |
|
|
476 ASSERT (col >= 0 && col < n_col) ; |
|
|
477 ASSERT (W [col] == Wflag) ; |
|
|
478 if (col > j2) |
|
|
479 { |
|
|
480 /* keep this col in the pattern of the new row j */ |
|
|
481 Ai [pdest++] = col ; |
|
|
482 if (col < parent) |
|
|
483 { |
|
|
484 parent = col ; |
|
|
485 } |
|
|
486 } |
|
|
487 } |
|
|
488 } |
|
|
489 |
|
|
490 DEBUG1 (("Parent ["ID"] = "ID"\n", j, parent)) ; |
|
|
491 ASSERT (parent > j2) ; |
|
|
492 |
|
|
493 if (parent == n_col) |
|
|
494 { |
|
|
495 /* this front has no parent - it is the root of a subtree */ |
|
|
496 parent = EMPTY ; |
|
|
497 } |
|
|
498 |
|
|
499 #ifndef NDEBUG |
|
|
500 DEBUG3 (("\nFinal row "ID" of U after supercolumn detection:\n", j)) ; |
|
|
501 for (p = Up [j] ; p < pdest ; p++) |
|
|
502 { |
|
|
503 ASSERT (Ai [p] >= 0 && Ai [p] < n_col) ; |
|
|
504 DEBUG3 ((" "ID" ("ID")", Ai [p], W [Ai [p]])) ; |
|
|
505 ASSERT (W [Ai [p]] == Wflag) ; |
|
|
506 } |
|
|
507 DEBUG3 (("\n")) ; |
|
|
508 #endif |
|
|
509 |
|
|
510 /* ================================================================== */ |
|
|
511 /* === frontal matrix =============================================== */ |
|
|
512 /* ================================================================== */ |
|
|
513 |
|
|
514 /* front has Front_npivcol [j] pivot columns */ |
|
|
515 /* entire front is Front_nrows [j] -by- Front_ncols [j] */ |
|
|
516 /* j is first column in the front */ |
|
|
517 |
|
|
518 npivcol = npivots ; |
|
|
519 fallrows = thickness ; |
|
|
520 fallcols = npivots + pdest - Up [j] ; |
|
|
521 |
|
|
522 /* number of pivots in the front (rows and columns) */ |
|
|
523 fpiv = MIN (npivcol, fallrows) ; |
|
|
524 |
|
|
525 /* size of contribution block */ |
|
|
526 frows = fallrows - fpiv ; |
|
|
527 fcols = fallcols - fpiv ; |
|
|
528 |
|
|
529 if (frows == 0 || fcols == 0) |
|
|
530 { |
|
|
531 /* front has no contribution block and thus needs no parent */ |
|
|
532 DEBUG1 (("Frontal matrix evaporation\n")) ; |
|
|
533 Up [j] = EMPTY ; |
|
|
534 parent = EMPTY ; |
|
|
535 } |
|
|
536 |
|
|
537 Front_npivcol [j] = npivots ; |
|
|
538 Front_nrows [j] = fallrows ; |
|
|
539 Front_ncols [j] = fallcols ; |
|
|
540 Front_parent [j] = parent ; |
|
|
541 ASSERT (npivots > 0) ; |
|
|
542 |
|
|
543 /* Front_parent [j] is the first column of the parent frontal matrix */ |
|
|
544 |
|
|
545 DEBUG1 (("\n\n==== Front "ID", nfr "ID" pivot columns "ID":"ID |
|
|
546 " all front: "ID"-by-"ID" Parent: "ID"\n", j, nfr, j,j+npivots-1, |
|
|
547 Front_nrows [j], Front_ncols [j], Front_parent [j])) ; |
|
|
548 nfr++ ; |
|
|
549 |
|
|
550 /* ================================================================== */ |
|
|
551 /* === prepare this row for its parent ============================== */ |
|
|
552 /* ================================================================== */ |
|
|
553 |
|
|
554 if (parent != EMPTY) |
|
|
555 { |
|
|
556 Link [j] = Link [parent] ; |
|
|
557 Link [parent] = j ; |
|
|
558 } |
|
|
559 |
|
|
560 ASSERT (jnext > j) ; |
|
|
561 |
|
|
562 jlast = j ; |
|
|
563 } |
|
|
564 |
|
|
565 /* ====================================================================== */ |
|
|
566 /* === postorder the fronts ============================================= */ |
|
|
567 /* ====================================================================== */ |
|
|
568 |
|
|
569 *nfr_out = nfr ; |
|
|
570 |
|
|
571 Front_order = W ; /* use W for Front_order [ */ |
|
|
572 |
|
|
573 if (fixQ) |
|
|
574 { |
|
|
575 /* do not postorder the fronts if Q is fixed */ |
|
|
576 DEBUG1 (("\nNo postorder (Q is fixed)\n")) ; |
|
|
577 k = 0 ; |
|
|
578 /* Pragma added May 14, 2003. The Intel compiler icl 6.0 (an old |
|
|
579 * version) incorrectly vectorizes this loop. */ |
|
|
580 #pragma novector |
|
|
581 for (j = 0 ; j < n_col ; j++) |
|
|
582 { |
|
|
583 if (Front_npivcol [j] > 0) |
|
|
584 { |
|
|
585 Front_order [j] = k++ ; |
|
|
586 DEBUG1 (("Front order of j: "ID" is:"ID"\n", j, |
|
|
587 Front_order [j])) ; |
|
|
588 } |
|
|
589 else |
|
|
590 { |
|
|
591 Front_order [j] = EMPTY ; |
|
|
592 } |
|
|
593 } |
|
|
594 } |
|
|
595 else |
|
|
596 { |
|
|
597 |
|
|
598 /* use Ap for Front_child and use Link for Front_sibling [ */ |
|
|
599 Front_child = Ap ; |
|
|
600 Front_sibling = Link ; |
|
|
601 |
|
|
602 /* use Ai for Front_stack, size of Ai is >= 2*n_col */ |
|
|
603 Front_stack = Ai ; |
|
|
604 Front_size = Front_stack + n_col ; |
|
|
605 |
|
|
606 UMF_fsize (n_col, Front_size, Front_nrows, Front_ncols, |
|
|
607 Front_parent, Front_npivcol) ; |
|
|
608 |
|
|
609 AMD_postorder (n_col, Front_parent, Front_npivcol, Front_size, |
|
|
610 Front_order, Front_child, Front_sibling, Front_stack) ; |
|
|
611 |
|
|
612 /* done with Front_child, Front_sibling, Front_size, and Front_stack ]*/ |
|
|
613 |
|
|
614 /* ------------------------------------------------------------------ */ |
|
|
615 /* construct the column permutation (return in Up) */ |
|
|
616 /* ------------------------------------------------------------------ */ |
|
|
617 |
|
|
618 /* Front_order [i] = k means that front i is kth front in the new order. |
|
|
619 * i is in the range 0 to n_col-1, and k is in the range 0 to nfr-1 */ |
|
|
620 |
|
|
621 /* Use Ai as workspace for Winv [ */ |
|
|
622 Winv = Ai ; |
|
|
623 for (k = 0 ; k < nfr ; k++) |
|
|
624 { |
|
|
625 Winv [k] = EMPTY ; |
|
|
626 } |
|
|
627 |
|
|
628 /* compute the inverse of Front_order, so that Winv [k] = i */ |
|
|
629 /* if Front_order [i] = k */ |
|
|
630 |
|
|
631 DEBUG1 (("\n\nComputing output column permutation:\n")) ; |
|
|
632 for (i = 0 ; i < n_col ; i++) |
|
|
633 { |
|
|
634 k = Front_order [i] ; |
|
|
635 if (k != EMPTY) |
|
|
636 { |
|
|
637 DEBUG1 (("Front "ID" new order: "ID"\n", i, k)) ; |
|
|
638 ASSERT (k >= 0 && k < nfr) ; |
|
|
639 ASSERT (Winv [k] == EMPTY) ; |
|
|
640 Winv [k] = i ; |
|
|
641 } |
|
|
642 } |
|
|
643 |
|
|
644 /* Use Up as output permutation */ |
|
|
645 kk = 0 ; |
|
|
646 for (k = 0 ; k < nfr ; k++) |
|
|
647 { |
|
|
648 i = Winv [k] ; |
|
|
649 DEBUG1 (("Old Front "ID" New Front "ID" npivots "ID" nrows "ID |
|
|
650 " ncols "ID"\n", |
|
|
651 i, k, Front_npivcol [i], Front_nrows [i], Front_ncols [i])) ; |
|
|
652 ASSERT (i >= 0 && i < n_col) ; |
|
|
653 ASSERT (Front_npivcol [i] > 0) ; |
|
|
654 for (npiv = 0 ; npiv < Front_npivcol [i] ; npiv++) |
|
|
655 { |
|
|
656 Up [kk] = i + npiv ; |
|
|
657 DEBUG1 ((" Cperm ["ID"] = "ID"\n", kk, Up [kk])) ; |
|
|
658 kk++ ; |
|
|
659 } |
|
|
660 } |
|
|
661 ASSERT (kk == n_col) ; |
|
|
662 |
|
|
663 /* Winv no longer needed ] */ |
|
|
664 } |
|
|
665 |
|
|
666 /* ---------------------------------------------------------------------- */ |
|
|
667 /* apply the postorder traversal to renumber the frontal matrices */ |
|
|
668 /* (or pack them in same order, if fixQ) */ |
|
|
669 /* ---------------------------------------------------------------------- */ |
|
|
670 |
|
|
671 /* use Ai as workspace */ |
|
|
672 |
|
|
673 UMF_apply_order (Front_npivcol, Front_order, Ai, n_col, nfr) ; |
|
|
674 UMF_apply_order (Front_nrows, Front_order, Ai, n_col, nfr) ; |
|
|
675 UMF_apply_order (Front_ncols, Front_order, Ai, n_col, nfr) ; |
|
|
676 UMF_apply_order (Front_parent, Front_order, Ai, n_col, nfr) ; |
|
|
677 |
|
|
678 /* fix the parent to refer to the new numbering */ |
|
|
679 for (i = 0 ; i < nfr ; i++) |
|
|
680 { |
|
|
681 parent = Front_parent [i] ; |
|
|
682 if (parent != EMPTY) |
|
|
683 { |
|
|
684 ASSERT (parent >= 0 && parent < n_col) ; |
|
|
685 ASSERT (Front_order [parent] >= 0 && Front_order [parent] < nfr) ; |
|
|
686 Front_parent [i] = Front_order [parent] ; |
|
|
687 } |
|
|
688 } |
|
|
689 |
|
|
690 /* Front_order longer needed ] */ |
|
|
691 |
|
|
692 #ifndef NDEBUG |
|
|
693 DEBUG1 (("\nFinal frontal matrices:\n")) ; |
|
|
694 for (i = 0 ; i < nfr ; i++) |
|
|
695 { |
|
|
696 DEBUG1 (("Final front "ID": npiv "ID" nrows "ID" ncols "ID" parent " |
|
|
697 ID"\n", i, Front_npivcol [i], Front_nrows [i], |
|
|
698 Front_ncols [i], Front_parent [i])) ; |
|
|
699 } |
|
|
700 #endif |
|
|
701 |
|
|
702 *p_ncompactions = ncompactions ; |
|
|
703 return (TRUE) ; |
|
|
704 } |
