1 /* ========================================================================== */
2 /* === UMF_singletons ======================================================= */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /* Find and order the row and column singletons of a matrix A. If there are
11 * row and column singletons, the output is a row and column permutation such
12 * that the matrix is in the following form:
13 *
14 * x x x x x x x x x
15 * . x x x x x x x x
16 * . . x x x x x x x
17 * . . . x . . . . .
18 * . . . x x . . . .
19 * . . . x x s s s s
20 * . . . x x s s s s
21 * . . . x x s s s s
22 * . . . x x s s s s
23 *
24 * The above example has 3 column singletons (the first three columns and
25 * their corresponding pivot rows) and 2 row singletons. The singletons are
26 * ordered first, because they have zero Markowitz cost. The LU factorization
27 * for these first five rows and columns is free - there is no work to do
28 * (except to scale the pivot columns for the 2 row singletons), and no
29 * fill-in occurs. The remaining submatrix (4-by-4 in the above example)
30 * has no rows or columns with degree one. It may have empty rows or columns.
31 *
32 * This algorithm does not perform a full permutation to block triangular
33 * form. If there are one or more singletons, then the matrix can be
34 * permuted to block triangular form, but UMFPACK does not perform the full
35 * BTF permutation (see also "dmperm" in MATLAB, CSparse cs_dmperm,
36 * and SuiteSparse/BTF).
37 */
38
39 #include "umf_internal.h"
40 #include "umf_singletons.h"
41
42 #ifndef NDEBUG
43
44 /* ========================================================================== */
45 /* === debug routines ======================================================= */
46 /* ========================================================================== */
47
48 /* Dump the singleton queue */
49
dump_singletons(Int head,Int tail,Int Next[],char * name,Int Deg[],Int n)50 PRIVATE void dump_singletons
51 (
52 Int head, /* head of the queue */
53 Int tail, /* tail of the queue */
54 Int Next [ ], /* Next [i] is the next object after i */
55 char *name, /* "row" or "col" */
56 Int Deg [ ], /* Deg [i] is the degree of object i */
57 Int n /* objects are in the range 0 to n-1 */
58 )
59 {
60 Int i, next, cnt ;
61 DEBUG6 (("%s Singleton list: head "ID" tail "ID"\n", name, head, tail)) ;
62 i = head ;
63 ASSERT (head >= EMPTY && head < n) ;
64 ASSERT (tail >= EMPTY && tail < n) ;
65 cnt = 0 ;
66 while (i != EMPTY)
67 {
68 DEBUG7 ((" "ID": "ID" deg: "ID"\n", cnt, i, Deg [i])) ;
69 ASSERT (i >= 0 && i < n) ;
70 next = Next [i] ;
71 if (i == tail) ASSERT (next == EMPTY) ;
72 i = next ;
73 cnt++ ;
74 ASSERT (cnt <= n) ;
75 }
76 }
77
dump_mat(char * xname,char * yname,Int nx,Int ny,const Int Xp[],const Int Xi[],Int Xdeg[],Int Ydeg[])78 PRIVATE void dump_mat
79 (
80 char *xname,
81 char *yname,
82 Int nx,
83 Int ny,
84 const Int Xp [ ],
85 const Int Xi [ ],
86 Int Xdeg [ ],
87 Int Ydeg [ ]
88 )
89 {
90 Int x, y, p, p1, p2, xdeg, do_xdeg, ydeg ;
91 DEBUG6 (("\n ==== Dump %s mat:\n", xname)) ;
92 for (x = 0 ; x < nx ; x++)
93 {
94 p1 = Xp [x] ;
95 p2 = Xp [x+1] ;
96 xdeg = Xdeg [x] ;
97 DEBUG6 (("Dump %s "ID" p1 "ID" p2 "ID" deg "ID"\n",
98 xname, x, p1, p2, xdeg)) ;
99 do_xdeg = (xdeg >= 0) ;
100 for (p = p1 ; p < p2 ; p++)
101 {
102 y = Xi [p] ;
103 DEBUG7 ((" %s "ID" deg: ", yname, y)) ;
104 ASSERT (y >= 0 && y < ny) ;
105 ydeg = Ydeg [y] ;
106 DEBUG7 ((ID"\n", ydeg)) ;
107 if (do_xdeg && ydeg >= 0)
108 {
109 xdeg-- ;
110 }
111 }
112 ASSERT (IMPLIES (do_xdeg, xdeg == 0)) ;
113 }
114 }
115 #endif
116
117 /* ========================================================================== */
118 /* === create_row_form ====================================================== */
119 /* ========================================================================== */
120
121 /* Create the row-form R of the column-form input matrix A. This could be done
122 * by UMF_transpose, except that Rdeg has already been computed.
123 */
124
create_row_form(Int n_row,Int n_col,const Int Ap[],const Int Ai[],Int Rdeg[],Int Rp[],Int Ri[],Int W[])125 PRIVATE void create_row_form
126 (
127 /* input, not modified: */
128 Int n_row, /* A is n_row-by-n_col, nz = Ap [n_col] */
129 Int n_col,
130 const Int Ap [ ], /* Ap [0..n_col]: column pointers for A */
131 const Int Ai [ ], /* Ai [0..nz-1]: row indices for A */
132 Int Rdeg [ ], /* Rdeg [0..n_row-1]: row degrees */
133
134 /* output, not defined on input: */
135 Int Rp [ ], /* Rp [0..n_row]: row pointers for R */
136 Int Ri [ ], /* Ri [0..nz-1]: column indices for R */
137
138 /* workspace, not defined on input or output */
139 Int W [ ] /* size n_row */
140 )
141 {
142 Int row, col, p, p2 ;
143
144 /* create the row pointers */
145 Rp [0] = 0 ;
146 W [0] = 0 ;
147 for (row = 0 ; row < n_row ; row++)
148 {
149 Rp [row+1] = Rp [row] + Rdeg [row] ;
150 W [row] = Rp [row] ;
151 }
152
153 /* create the indices for the row-form */
154 for (col = 0 ; col < n_col ; col++)
155 {
156 p2 = Ap [col+1] ;
157 for (p = Ap [col] ; p < p2 ; p++)
158 {
159 Ri [W [Ai [p]]++] = col ;
160 }
161 }
162 }
163
164 /* ========================================================================== */
165 /* === order_singletons ===================================================== */
166 /* ========================================================================== */
167
order_singletons(Int k,Int head,Int tail,Int Next[],Int Xdeg[],Int Xperm[],const Int Xp[],const Int Xi[],Int Ydeg[],Int Yperm[],const Int Yp[],const Int Yi[],char * xname,char * yname,Int nx,Int ny)168 PRIVATE int order_singletons /* return new number of singletons */
169 (
170 Int k, /* the number of singletons so far */
171 Int head,
172 Int tail,
173 Int Next [ ],
174 Int Xdeg [ ], Int Xperm [ ], const Int Xp [ ], const Int Xi [ ],
175 Int Ydeg [ ], Int Yperm [ ], const Int Yp [ ], const Int Yi [ ]
176 #ifndef NDEBUG
177 , char *xname, char *yname, Int nx, Int ny
178 #endif
179 )
180 {
181 Int xpivot, x, y, ypivot, p, p2, deg ;
182
183 #ifndef NDEBUG
184 Int i, k1 = k ;
185 dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
186 dump_mat (xname, yname, nx, ny, Xp, Xi, Xdeg, Ydeg) ;
187 dump_mat (yname, xname, ny, nx, Yp, Yi, Ydeg, Xdeg) ;
188 #endif
189
190 while (head != EMPTY)
191 {
192 /* remove the singleton at the head of the queue */
193 xpivot = head ;
194 DEBUG1 (("------ Order %s singleton: "ID"\n", xname, xpivot)) ;
195 head = Next [xpivot] ;
196 if (head == EMPTY) tail = EMPTY ;
197
198 #ifndef NDEBUG
199 if (k % 100 == 0) dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
200 #endif
201
202 ASSERT (Xdeg [xpivot] >= 0) ;
203 if (Xdeg [xpivot] != 1)
204 {
205 /* This row/column x is empty. The matrix is singular.
206 * x will be ordered last in Xperm. */
207 DEBUG1 (("empty %s, after singletons removed\n", xname)) ;
208 continue ;
209 }
210
211 /* find the ypivot to match with this xpivot */
212 #ifndef NDEBUG
213 /* there can only be one ypivot, since the degree of x is 1 */
214 deg = 0 ;
215 p2 = Xp [xpivot+1] ;
216 for (p = Xp [xpivot] ; p < p2 ; p++)
217 {
218 y = Xi [p] ;
219 DEBUG1 (("%s: "ID"\n", yname, y)) ;
220 if (Ydeg [y] >= 0)
221 {
222 /* this is a live index in this xpivot vector */
223 deg++ ;
224 }
225 }
226 ASSERT (deg == 1) ;
227 #endif
228
229 ypivot = EMPTY ;
230 p2 = Xp [xpivot+1] ;
231 for (p = Xp [xpivot] ; p < p2 ; p++)
232 {
233 y = Xi [p] ;
234 DEBUG1 (("%s: "ID"\n", yname, y)) ;
235 if (Ydeg [y] >= 0)
236 {
237 /* this is a live index in this xpivot vector */
238 ypivot = y ;
239 break ;
240 }
241 }
242
243 DEBUG1 (("Pivot %s: "ID"\n", yname, ypivot)) ;
244 ASSERT (ypivot != EMPTY) ;
245 DEBUG1 (("deg "ID"\n", Ydeg [ypivot])) ;
246 ASSERT (Ydeg [ypivot] >= 0) ;
247
248 /* decrement the degrees after removing this singleton */
249 DEBUG1 (("p1 "ID"\n", Yp [ypivot])) ;
250 DEBUG1 (("p2 "ID"\n", Yp [ypivot+1])) ;
251 p2 = Yp [ypivot+1] ;
252 for (p = Yp [ypivot] ; p < p2 ; p++)
253 {
254 x = Yi [p] ;
255 DEBUG1 ((" %s: "ID" deg: "ID"\n", xname, x, Xdeg [x])) ;
256 if (Xdeg [x] < 0) continue ;
257 ASSERT (Xdeg [x] > 0) ;
258 if (x == xpivot) continue ;
259 deg = --(Xdeg [x]) ;
260 ASSERT (Xdeg [x] >= 0) ;
261 if (deg == 1)
262 {
263 /* this is a new singleton, put at the end of the queue */
264 Next [x] = EMPTY ;
265 if (head == EMPTY)
266 {
267 head = x ;
268 }
269 else
270 {
271 ASSERT (tail != EMPTY) ;
272 Next [tail] = x ;
273 }
274 tail = x ;
275 DEBUG1 ((" New %s singleton: "ID"\n", xname, x)) ;
276 #ifndef NDEBUG
277 if (k % 100 == 0)
278 {
279 dump_singletons (head, tail, Next, xname, Xdeg, nx) ;
280 }
281 #endif
282 }
283 }
284
285 /* flag the xpivot and ypivot by FLIP'ing the degrees */
286 Xdeg [xpivot] = FLIP (1) ;
287 Ydeg [ypivot] = FLIP (Ydeg [ypivot]) ;
288
289 /* keep track of the pivot row and column */
290 Xperm [k] = xpivot ;
291 Yperm [k] = ypivot ;
292 k++ ;
293
294 #ifndef NDEBUG
295 if (k % 1000 == 0)
296 {
297 dump_mat (xname, yname, nx, ny, Xp, Xi, Xdeg, Ydeg) ;
298 dump_mat (yname, xname, ny, nx, Yp, Yi, Ydeg, Xdeg) ;
299 }
300 #endif
301 }
302
303 #ifndef NDEBUG
304 DEBUGm4 (("%s singletons: k = "ID"\n", xname, k)) ;
305 for (i = k1 ; i < k ; i++)
306 {
307 DEBUG1 ((" %s: "ID" %s: "ID"\n", xname, Xperm [i], yname, Yperm [i])) ;
308 }
309 ASSERT (k > 0) ;
310 #endif
311
312 return (k) ;
313 }
314
315 /* ========================================================================== */
316 /* === find_any_singletons ================================================== */
317 /* ========================================================================== */
318
find_any_singletons(Int n_row,Int n_col,const Int Ap[],const Int Ai[],Int Cdeg[],Int Rdeg[],Int Cperm[],Int Rperm[],Int * p_n1r,Int * p_n1c,Int Rp[],Int Ri[],Int W[],Int Next[])319 PRIVATE Int find_any_singletons /* returns # of singletons found */
320 (
321 /* input, not modified: */
322 Int n_row,
323 Int n_col,
324 const Int Ap [ ], /* size n_col+1 */
325 const Int Ai [ ], /* size nz = Ap [n_col] */
326
327 /* input, modified on output: */
328 Int Cdeg [ ], /* size n_col */
329 Int Rdeg [ ], /* size n_row */
330
331 /* output, not defined on input: */
332 Int Cperm [ ], /* size n_col */
333 Int Rperm [ ], /* size n_row */
334 Int *p_n1r, /* # of row singletons */
335 Int *p_n1c, /* # of col singletons */
336
337 /* workspace, not defined on input or output */
338 Int Rp [ ], /* size n_row+1 */
339 Int Ri [ ], /* size nz */
340 Int W [ ], /* size n_row */
341 Int Next [ ] /* size MAX (n_row, n_col) */
342 )
343 {
344 Int n1, col, row, row_form, head, tail, n1r, n1c ;
345
346 /* ---------------------------------------------------------------------- */
347 /* eliminate column singletons */
348 /* ---------------------------------------------------------------------- */
349
350 n1 = 0 ;
351 n1r = 0 ;
352 n1c = 0 ;
353 row_form = FALSE ;
354
355 head = EMPTY ;
356 tail = EMPTY ;
357 for (col = n_col-1 ; col >= 0 ; col--)
358 {
359 if (Cdeg [col] == 1)
360 {
361 /* put the column singleton in the queue */
362 if (head == EMPTY) tail = col ;
363 Next [col] = head ;
364 head = col ;
365 DEBUG1 (("Column singleton: "ID"\n", col)) ;
366 }
367 }
368
369 if (head != EMPTY)
370 {
371
372 /* ------------------------------------------------------------------ */
373 /* create the row-form of A */
374 /* ------------------------------------------------------------------ */
375
376 create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
377 row_form = TRUE ;
378
379 /* ------------------------------------------------------------------ */
380 /* find and order the column singletons */
381 /* ------------------------------------------------------------------ */
382
383 n1 = order_singletons (0, head, tail, Next,
384 Cdeg, Cperm, Ap, Ai,
385 Rdeg, Rperm, Rp, Ri
386 #ifndef NDEBUG
387 , "col", "row", n_col, n_row
388 #endif
389 ) ;
390 n1c = n1 ;
391 }
392
393 /* ---------------------------------------------------------------------- */
394 /* eliminate row singletons */
395 /* ---------------------------------------------------------------------- */
396
397 head = EMPTY ;
398 tail = EMPTY ;
399 for (row = n_row-1 ; row >= 0 ; row--)
400 {
401 if (Rdeg [row] == 1)
402 {
403 /* put the row singleton in the queue */
404 if (head == EMPTY) tail = row ;
405 Next [row] = head ;
406 head = row ;
407 DEBUG1 (("Row singleton: "ID"\n", row)) ;
408 }
409 }
410
411 if (head != EMPTY)
412 {
413
414 /* ------------------------------------------------------------------ */
415 /* create the row-form of A, if not already created */
416 /* ------------------------------------------------------------------ */
417
418 if (!row_form)
419 {
420 create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
421 }
422
423 /* ------------------------------------------------------------------ */
424 /* find and order the row singletons */
425 /* ------------------------------------------------------------------ */
426
427 n1 = order_singletons (n1, head, tail, Next,
428 Rdeg, Rperm, Rp, Ri,
429 Cdeg, Cperm, Ap, Ai
430 #ifndef NDEBUG
431 , "row", "col", n_row, n_col
432 #endif
433 ) ;
434 n1r = n1 - n1c ;
435 }
436
437 DEBUG0 (("n1 "ID"\n", n1)) ;
438 *p_n1r = n1r ;
439 *p_n1c = n1c ;
440 return (n1) ;
441 }
442
443 /* ========================================================================== */
444 /* === find_user_singletons ================================================= */
445 /* ========================================================================== */
446
find_user_singletons(Int n_row,Int n_col,const Int Ap[],const Int Ai[],const Int Quser[],Int Cdeg[],Int Rdeg[],Int Cperm[],Int Rperm[],Int * p_n1r,Int * p_n1c,Int Rp[],Int Ri[],Int W[])447 PRIVATE Int find_user_singletons /* returns # singletons found */
448 (
449 /* input, not modified: */
450 Int n_row,
451 Int n_col,
452 const Int Ap [ ], /* size n_col+1 */
453 const Int Ai [ ], /* size nz = Ap [n_col] */
454 const Int Quser [ ], /* size n_col if present */
455
456 /* input, modified on output: */
457 Int Cdeg [ ], /* size n_col */
458 Int Rdeg [ ], /* size n_row */
459
460 /* output, not defined on input */
461 Int Cperm [ ], /* size n_col */
462 Int Rperm [ ], /* size n_row */
463 Int *p_n1r, /* # of row singletons */
464 Int *p_n1c, /* # of col singletons */
465
466 /* workspace, not defined on input or output */
467 Int Rp [ ], /* size n_row+1 */
468 Int Ri [ ], /* size nz */
469 Int W [ ] /* size n_row */
470 )
471 {
472 Int n1, col, row, p, p2, pivcol, pivrow, found, k, n1r, n1c ;
473
474 n1 = 0 ;
475 n1r = 0 ;
476 n1c = 0 ;
477 *p_n1r = 0 ;
478 *p_n1c = 0 ;
479
480 /* find singletons in the user column permutation, Quser */
481 pivcol = Quser [0] ;
482 found = (Cdeg [pivcol] == 1) ;
483 DEBUG0 (("Is first col: "ID" a col singleton?: "ID"\n", pivcol, found)) ;
484 if (!found)
485 {
486 /* the first column is not a column singleton, check for a row
487 * singleton in the first column. */
488 for (p = Ap [pivcol] ; p < Ap [pivcol+1] ; p++)
489 {
490 if (Rdeg [Ai [p]] == 1)
491 {
492 DEBUG0 (("Row singleton in first col: "ID" row: "ID"\n",
493 pivcol, Ai [p])) ;
494 found = TRUE ;
495 break ;
496 }
497 }
498 }
499
500 if (!found)
501 {
502 /* no singletons in the leading part of A (:,Quser) */
503 return (0) ;
504 }
505
506 /* there is at least one row or column singleton. Look for more. */
507 create_row_form (n_row, n_col, Ap, Ai, Rdeg, Rp, Ri, W) ;
508
509 n1 = 0 ;
510
511 for (k = 0 ; k < n_col ; k++)
512 {
513 pivcol = Quser [k] ;
514 pivrow = EMPTY ;
515
516 /* ------------------------------------------------------------------ */
517 /* check if col is a column singleton, or contains a row singleton */
518 /* ------------------------------------------------------------------ */
519
520 found = (Cdeg [pivcol] == 1) ;
521
522 if (found)
523 {
524
525 /* -------------------------------------------------------------- */
526 /* pivcol is a column singleton */
527 /* -------------------------------------------------------------- */
528
529 DEBUG0 (("Found a col singleton: k "ID" pivcol "ID"\n", k, pivcol));
530
531 /* find the pivrow to match with this pivcol */
532 #ifndef NDEBUG
533 /* there can only be one pivrow, since the degree of pivcol is 1 */
534 {
535 Int deg = 0 ;
536 p2 = Ap [pivcol+1] ;
537 for (p = Ap [pivcol] ; p < p2 ; p++)
538 {
539 row = Ai [p] ;
540 DEBUG1 (("row: "ID"\n", row)) ;
541 if (Rdeg [row] >= 0)
542 {
543 /* this is a live index in this column vector */
544 deg++ ;
545 }
546 }
547 ASSERT (deg == 1) ;
548 }
549 #endif
550
551 p2 = Ap [pivcol+1] ;
552 for (p = Ap [pivcol] ; p < p2 ; p++)
553 {
554 row = Ai [p] ;
555 DEBUG1 (("row: "ID"\n", row)) ;
556 if (Rdeg [row] >= 0)
557 {
558 /* this is a live index in this pivcol vector */
559 pivrow = row ;
560 break ;
561 }
562 }
563
564 DEBUG1 (("Pivot row: "ID"\n", pivrow)) ;
565 ASSERT (pivrow != EMPTY) ;
566 DEBUG1 (("deg "ID"\n", Rdeg [pivrow])) ;
567 ASSERT (Rdeg [pivrow] >= 0) ;
568
569 /* decrement the degrees after removing this col singleton */
570 DEBUG1 (("p1 "ID"\n", Rp [pivrow])) ;
571 DEBUG1 (("p2 "ID"\n", Rp [pivrow+1])) ;
572 p2 = Rp [pivrow+1] ;
573 for (p = Rp [pivrow] ; p < p2 ; p++)
574 {
575 col = Ri [p] ;
576 DEBUG1 ((" col: "ID" deg: "ID"\n", col, Cdeg [col])) ;
577 if (Cdeg [col] < 0) continue ;
578 ASSERT (Cdeg [col] > 0) ;
579 Cdeg [col]-- ;
580 ASSERT (Cdeg [col] >= 0) ;
581 }
582
583 /* flag the pivcol and pivrow by FLIP'ing the degrees */
584 Cdeg [pivcol] = FLIP (1) ;
585 Rdeg [pivrow] = FLIP (Rdeg [pivrow]) ;
586 n1c++ ;
587
588 }
589 else
590 {
591
592 /* -------------------------------------------------------------- */
593 /* pivcol may contain a row singleton */
594 /* -------------------------------------------------------------- */
595
596 p2 = Ap [pivcol+1] ;
597 for (p = Ap [pivcol] ; p < p2 ; p++)
598 {
599 pivrow = Ai [p] ;
600 if (Rdeg [pivrow] == 1)
601 {
602 DEBUG0 (("Row singleton in pivcol: "ID" row: "ID"\n",
603 pivcol, pivrow)) ;
604 found = TRUE ;
605 break ;
606 }
607 }
608
609 if (!found)
610 {
611 DEBUG0 (("End of user singletons\n")) ;
612 break ;
613 }
614
615 #ifndef NDEBUG
616 /* there can only be one pivrow, since the degree of pivcol is 1 */
617 {
618 Int deg = 0 ;
619 p2 = Rp [pivrow+1] ;
620 for (p = Rp [pivrow] ; p < p2 ; p++)
621 {
622 col = Ri [p] ;
623 DEBUG1 (("col: "ID" cdeg::: "ID"\n", col, Cdeg [col])) ;
624 if (Cdeg [col] >= 0)
625 {
626 /* this is a live index in this column vector */
627 ASSERT (col == pivcol) ;
628 deg++ ;
629 }
630 }
631 ASSERT (deg == 1) ;
632 }
633 #endif
634
635 DEBUG1 (("Pivot row: "ID"\n", pivrow)) ;
636 DEBUG1 (("pivcol deg "ID"\n", Cdeg [pivcol])) ;
637 ASSERT (Cdeg [pivcol] > 1) ;
638
639 /* decrement the degrees after removing this row singleton */
640 DEBUG1 (("p1 "ID"\n", Ap [pivcol])) ;
641 DEBUG1 (("p2 "ID"\n", Ap [pivcol+1])) ;
642 p2 = Ap [pivcol+1] ;
643 for (p = Ap [pivcol] ; p < p2 ; p++)
644 {
645 row = Ai [p] ;
646 DEBUG1 ((" row: "ID" deg: "ID"\n", row, Rdeg [row])) ;
647 if (Rdeg [row] < 0) continue ;
648 ASSERT (Rdeg [row] > 0) ;
649 Rdeg [row]-- ;
650 ASSERT (Rdeg [row] >= 0) ;
651 }
652
653 /* flag the pivcol and pivrow by FLIP'ing the degrees */
654 Cdeg [pivcol] = FLIP (Cdeg [pivcol]) ;
655 Rdeg [pivrow] = FLIP (1) ;
656 n1r++ ;
657 }
658
659 /* keep track of the pivot row and column */
660 Cperm [k] = pivcol ;
661 Rperm [k] = pivrow ;
662 n1++ ;
663
664 #ifndef NDEBUG
665 dump_mat ("col", "row", n_col, n_row, Ap, Ai, Cdeg, Rdeg) ;
666 dump_mat ("row", "col", n_row, n_col, Rp, Ri, Rdeg, Cdeg) ;
667 #endif
668
669 }
670
671 DEBUGm4 (("User singletons found: "ID"\n", n1)) ;
672 ASSERT (n1 > 0) ;
673
674 *p_n1r = n1r ;
675 *p_n1c = n1c ;
676 return (n1) ;
677 }
678
679 /* ========================================================================== */
680 /* === finish_permutation =================================================== */
681 /* ========================================================================== */
682
683 /* Complete the permutation for the pruned submatrix. The singletons are
684 * already ordered, but remove their flags. Place rows/columns that are empty
685 * in the pruned submatrix at the end of the output permutation. This can only
686 * occur if the matrix is singular.
687 */
688
finish_permutation(Int n1,Int nx,Int Xdeg[],const Int Xuser[],Int Xperm[],Int * p_max_deg)689 PRIVATE Int finish_permutation
690 (
691 Int n1,
692 Int nx,
693 Int Xdeg [ ],
694 const Int Xuser [ ],
695 Int Xperm [ ],
696 Int *p_max_deg
697 )
698 {
699 Int nempty, x, deg, s, max_deg, k ;
700 nempty = 0 ;
701 s = n1 ;
702 max_deg = 0 ;
703 DEBUG0 (("n1 "ID" nempty "ID"\n", n1, nempty)) ;
704 for (k = 0 ; k < nx ; k++)
705 {
706 x = (Xuser != (Int *) NULL) ? Xuser [k] : k ;
707 DEBUG0 (("finish perm k "ID" x "ID" nx "ID"\n", k, x, nx)) ;
708 deg = Xdeg [x] ;
709 if (deg == 0)
710 {
711 /* this row/col is empty in the pruned submatrix */
712 ASSERT (s < nx - nempty) ;
713 DEBUG0 (("empty k "ID"\n", k)) ;
714 nempty++ ;
715 Xperm [nx - nempty] = x ;
716 }
717 else if (deg > 0)
718 {
719 /* this row/col is nonempty in the pruned submatrix */
720 ASSERT (s < nx - nempty) ;
721 Xperm [s++] = x ;
722 max_deg = MAX (max_deg, deg) ;
723 }
724 else
725 {
726 /* This is a singleton row/column - it is already ordered.
727 * Just clear the flag. */
728 Xdeg [x] = FLIP (deg) ;
729 }
730 }
731 ASSERT (s == nx - nempty) ;
732 *p_max_deg = max_deg ;
733 return (nempty) ;
734 }
735
736 /* ========================================================================== */
737 /* === UMF_singletons ======================================================= */
738 /* ========================================================================== */
739
UMF_singletons(Int n_row,Int n_col,const Int Ap[],const Int Ai[],const Int Quser[],Int strategy,Int do_singletons,Int Cdeg[],Int Cperm[],Int Rdeg[],Int Rperm[],Int InvRperm[],Int * p_n1,Int * p_n1c,Int * p_n1r,Int * p_nempty_col,Int * p_nempty_row,Int * p_is_sym,Int * p_max_rdeg,Int Rp[],Int Ri[],Int W[],Int Next[])740 GLOBAL Int UMF_singletons
741 (
742
743 /* input, not modified: */
744 Int n_row,
745 Int n_col,
746 const Int Ap [ ], /* size n_col+1 */
747 const Int Ai [ ], /* size nz = Ap [n_col] */
748 const Int Quser [ ], /* size n_col if present */
749 Int strategy, /* strategy requested by user */
750 Int do_singletons, /* if false, then do not look for singletons */
751
752 /* output, not defined on input: */
753 Int Cdeg [ ], /* size n_col */
754 Int Cperm [ ], /* size n_col */
755 Int Rdeg [ ], /* size n_row */
756 Int Rperm [ ], /* size n_row */
757 Int InvRperm [ ], /* size n_row, the inverse of Rperm */
758 Int *p_n1, /* # of col and row singletons */
759 Int *p_n1c, /* # of col singletons */
760 Int *p_n1r, /* # of row singletons */
761 Int *p_nempty_col, /* # of empty columns in pruned submatrix */
762 Int *p_nempty_row, /* # of empty columns in pruned submatrix */
763 Int *p_is_sym, /* TRUE if pruned submatrix is square and has been
764 * symmetrically permuted by Cperm and Rperm */
765 Int *p_max_rdeg, /* maximum Rdeg in pruned submatrix */
766
767 /* workspace, not defined on input or output */
768 Int Rp [ ], /* size n_row+1 */
769 Int Ri [ ], /* size nz */
770 Int W [ ], /* size n_row */
771 Int Next [ ] /* size MAX (n_row, n_col) */
772 )
773 {
774 Int n1, s, col, row, p, p1, p2, cdeg, last_row, is_sym, k,
775 nempty_row, nempty_col, max_cdeg, max_rdeg, n1c, n1r ;
776
777 /* ---------------------------------------------------------------------- */
778 /* initializations */
779 /* ---------------------------------------------------------------------- */
780
781 #ifndef NDEBUG
782 UMF_dump_start ( ) ;
783 DEBUGm4 (("Starting umf_singletons\n")) ;
784 #endif
785
786 /* ---------------------------------------------------------------------- */
787 /* scan the columns, check for errors and count row degrees */
788 /* ---------------------------------------------------------------------- */
789
790 if (Ap [0] != 0 || Ap [n_col] < 0)
791 {
792 return (UMFPACK_ERROR_invalid_matrix) ;
793 }
794 for (row = 0 ; row < n_row ; row++)
795 {
796 Rdeg [row] = 0 ;
797 }
798 for (col = 0 ; col < n_col ; col++)
799 {
800 p1 = Ap [col] ;
801 p2 = Ap [col+1] ;
802 cdeg = p2 - p1 ;
803 if (cdeg < 0)
804 {
805 return (UMFPACK_ERROR_invalid_matrix) ;
806 }
807 last_row = EMPTY ;
808 for (p = p1 ; p < p2 ; p++)
809 {
810 row = Ai [p] ;
811 if (row <= last_row || row >= n_row)
812 {
813 return (UMFPACK_ERROR_invalid_matrix) ;
814 }
815 Rdeg [row]++ ;
816 last_row = row ;
817 }
818 Cdeg [col] = cdeg ;
819 }
820
821 /* ---------------------------------------------------------------------- */
822 /* find singletons */
823 /* ---------------------------------------------------------------------- */
824
825 if (!do_singletons)
826 {
827 /* do not look for singletons at all */
828 n1 = 0 ;
829 n1r = 0 ;
830 n1c = 0 ;
831 }
832 else if (Quser != (Int *) NULL)
833 {
834 /* user has provided an input column ordering */
835 if (strategy == UMFPACK_STRATEGY_UNSYMMETRIC)
836 {
837 /* look for singletons, but respect the user's input permutation */
838 n1 = find_user_singletons (n_row, n_col, Ap, Ai, Quser,
839 Cdeg, Rdeg, Cperm, Rperm, &n1r, &n1c, Rp, Ri, W) ;
840 }
841 else
842 {
843 /* do not look for singletons if Quser given and strategy is
844 * not unsymmetric */
845 n1 = 0 ;
846 n1r = 0 ;
847 n1c = 0 ;
848 }
849 }
850 else
851 {
852 /* look for singletons anywhere */
853 n1 = find_any_singletons (n_row, n_col, Ap, Ai,
854 Cdeg, Rdeg, Cperm, Rperm, &n1r, &n1c, Rp, Ri, W, Next) ;
855 }
856
857 /* ---------------------------------------------------------------------- */
858 /* eliminate empty columns and complete the column permutation */
859 /* ---------------------------------------------------------------------- */
860
861 nempty_col = finish_permutation (n1, n_col, Cdeg, Quser, Cperm, &max_cdeg) ;
862
863 /* ---------------------------------------------------------------------- */
864 /* eliminate empty rows and complete the row permutation */
865 /* ---------------------------------------------------------------------- */
866
867 if (Quser != (Int *) NULL && strategy == UMFPACK_STRATEGY_SYMMETRIC)
868 {
869 /* rows should be symmetrically permuted according to Quser */
870 ASSERT (n_row == n_col) ;
871 nempty_row = finish_permutation (n1, n_row, Rdeg, Quser, Rperm,
872 &max_rdeg) ;
873 }
874 else
875 {
876 /* rows should not be symmetrically permuted according to Quser */
877 nempty_row = finish_permutation (n1, n_row, Rdeg, (Int *) NULL, Rperm,
878 &max_rdeg) ;
879 }
880
881 /* ---------------------------------------------------------------------- */
882 /* compute the inverse of Rperm */
883 /* ---------------------------------------------------------------------- */
884
885 for (k = 0 ; k < n_row ; k++)
886 {
887 ASSERT (Rperm [k] >= 0 && Rperm [k] < n_row) ;
888 InvRperm [Rperm [k]] = k ;
889 }
890
891 /* ---------------------------------------------------------------------- */
892 /* see if pruned submatrix is square and has been symmetrically permuted */
893 /* ---------------------------------------------------------------------- */
894
895 /* The prior version of this code (with a "break" statement; UMFPACK 5.2)
896 * causes UMFPACK to fail when optimization is enabled with gcc version
897 * 4.2.4 in a 64-bit Linux environment. The bug is a compiler bug, not a
898 * an UMFPACK bug. It is fixed in gcc version 4.3.2. However, as a
899 * workaround for the compiler, the code below has been "fixed". */
900
901 if (n_row == n_col && nempty_row == nempty_col)
902 {
903 /* is_sym is true if the submatrix is square, and
904 * Rperm [n1..n_row-nempty_row-1] = Cperm [n1..n_col-nempty_col-1] */
905 is_sym = TRUE ;
906 for (s = n1 ; /* replaced the break with this test: */ is_sym &&
907 /* the remainder of this test is unchanged from v5.2.0: */
908 s < n_col - nempty_col ; s++)
909 {
910 if (Cperm [s] != Rperm [s])
911 {
912 is_sym = FALSE ;
913 /* removed a break statement here, which is OK but it tickles
914 * the gcc 4.2.{3,4} compiler bug */
915 }
916 }
917 }
918 else
919 {
920 is_sym = FALSE ;
921 }
922
923 DEBUGm4 (("Submatrix square and symmetrically permuted? "ID"\n", is_sym)) ;
924 DEBUGm4 (("singletons "ID" row "ID" col "ID"\n", n1, n1r, n1c)) ;
925 DEBUGm4 (("Empty cols "ID" rows "ID"\n", nempty_col, nempty_row)) ;
926 *p_n1 = n1 ;
927 *p_n1r = n1r ;
928 *p_n1c = n1c ;
929 *p_is_sym = is_sym ;
930 *p_nempty_col = nempty_col ;
931 *p_nempty_row = nempty_row ;
932 *p_max_rdeg = max_rdeg ;
933 return (UMFPACK_OK) ;
934 }
935