1 /* ========================================================================== */
2 /* === UMF_colamd =========================================================== */
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 /*
11 UMF_colamd: an approximate minimum degree column ordering algorithm,
12 used as a preordering for UMFPACK.
13
14 NOTE: if this routine is used outside of UMFPACK, for a sparse Cholesky
15 factorization of (AQ)'*(AQ) or a QR factorization of A, then one line should
16 be removed (the "&& pivot_row_thickness > 0" expression). See the comment
17 regarding the Cholesky factorization, below.
18
19 Purpose:
20
21 Colamd computes a permutation Q such that the Cholesky factorization of
22 (AQ)'(AQ) has less fill-in and requires fewer floating point operations
23 than A'A. This also provides a good ordering for sparse partial
24 pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
25 factorization, and P is computed during numerical factorization via
26 conventional partial pivoting with row interchanges. Colamd is the
27 column ordering method used in SuperLU, part of the ScaLAPACK library.
28 It is also available as built-in function in MATLAB Version 6,
29 available from MathWorks, Inc. (http://www.mathworks.com). This
30 routine can be used in place of colmmd in MATLAB. By default, the \
31 and / operators in MATLAB perform a column ordering (using colmmd
32 or colamd) prior to LU factorization using sparse partial pivoting,
33 in the built-in MATLAB lu(A) routine.
34
35 This code is derived from Colamd Version 2.0.
36
37 Authors:
38
39 The authors of the COLAMD code itself are Stefan I. Larimore and Timothy A.
40 Davis. The algorithm was developed in collaboration
41 with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory.
42 The AMD metric on which this is based is by Patrick Amestoy, T. Davis,
43 and Iain Duff.
44
45 Date:
46
47 UMFPACK Version: see above.
48 COLAMD Version 2.0 was released on January 31, 2000.
49
50 Acknowledgements:
51
52 This work was supported by the National Science Foundation, under
53 grants DMS-9504974, DMS-9803599, and CCR-0203270.
54
55 UMFPACK: Copyright (c) 2003 by Timothy A. Davis. All Rights Reserved.
56
57 See the UMFPACK README file for the License for your use of this code.
58
59 Availability:
60
61 Both UMFPACK and the original unmodified colamd/symamd library are
62 available at http://www.suitesparse.com.
63
64 Changes for inclusion in UMFPACK:
65
66 * symamd, symamd_report, and colamd_report removed
67
68 * additional terms added to RowInfo, ColInfo, and stats
69
70 * Frontal matrix information computed for UMFPACK
71
72 * routines renamed
73
74 * column elimination tree post-ordering incorporated. In the original
75 version 2.0, this was performed in colamd.m.
76
77 For more information, see:
78
79 Amestoy, P. R. and Davis, T. A. and Duff, I. S.,
80 An approximate minimum degree ordering algorithm,
81 SIAM J. Matrix Analysis and Applic, vol 17, no 4., pp 886-905, 1996.
82
83 Davis, T. A. and Gilbert, J. R. and Larimore, S. I. and Ng, E. G.,
84 A column approximate minimum degree ordering algorithm,
85 ACM Trans. Math. Softw., vol 3, no 3, 2004
86
87 */
88
89 /* ========================================================================== */
90 /* === Description of user-callable routines ================================ */
91 /* ========================================================================== */
92
93 /*
94 ----------------------------------------------------------------------------
95 colamd_recommended: removed for UMFPACK
96 ----------------------------------------------------------------------------
97
98 ----------------------------------------------------------------------------
99 colamd_set_defaults:
100 ----------------------------------------------------------------------------
101
102 C syntax:
103
104 #include "colamd.h"
105 colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
106
107 Purpose:
108
109 Sets the default parameters. The use of this routine is optional.
110
111 Arguments:
112
113 double knobs [COLAMD_KNOBS] ; Output only.
114
115 Let c = knobs [COLAMD_DENSE_COL], r = knobs [COLAMD_DENSE_ROW].
116 Colamd: rows with more than max (16, r*16*sqrt(n_col))
117 entries are removed prior to ordering. Columns with more than
118 max (16, c*16*sqrt(n_row)) entries are removed prior to
119 ordering, and placed last in the output column ordering.
120
121 Symamd: removed for UMFPACK.
122
123 COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
124 respectively, in colamd.h. Default values of these two knobs
125 are both 0.5. Currently, only knobs [0] and knobs [1] are
126 used, but future versions may use more knobs. If so, they will
127 be properly set to their defaults by the future version of
128 colamd_set_defaults, so that the code that calls colamd will
129 not need to change, assuming that you either use
130 colamd_set_defaults, or pass a (double *) NULL pointer as the
131 knobs array to colamd or symamd.
132
133 knobs [COLAMD_AGGRESSIVE]: if nonzero, then perform aggressive
134 absorption. Otherwise, do not. This version does aggressive
135 absorption by default. COLAMD v2.1 (in MATLAB) always
136 does aggressive absorption (it doesn't have an option to turn
137 it off).
138
139 ----------------------------------------------------------------------------
140 colamd:
141 ----------------------------------------------------------------------------
142
143 C syntax:
144
145 #include "colamd.h"
146 Int UMF_colamd (Int n_row, Int n_col, Int Alen, Int *A, Int *p,
147 double knobs [COLAMD_KNOBS], Int stats [COLAMD_STATS]) ;
148
149 Purpose:
150
151 Computes a column ordering (Q) of A such that P(AQ)=LU or
152 (AQ)'AQ=LL' have less fill-in and require fewer floating point
153 operations than factorizing the unpermuted matrix A or A'A,
154 respectively.
155
156 Returns:
157
158 TRUE (1) if successful, FALSE (0) otherwise.
159
160 Arguments:
161
162 Int n_row ; Input argument.
163
164 Number of rows in the matrix A.
165 Restriction: n_row >= 0.
166 Colamd returns FALSE if n_row is negative.
167
168 Int n_col ; Input argument.
169
170 Number of columns in the matrix A.
171 Restriction: n_col >= 0.
172 Colamd returns FALSE if n_col is negative.
173
174 Int Alen ; Input argument.
175
176 Restriction (see note):
177 Alen >= 2*nnz + 8*(n_col+1) + 6*(n_row+1) + n_col
178 Colamd returns FALSE if these conditions are not met.
179
180 Note: this restriction makes an modest assumption regarding
181 the size of the two typedef's structures in colamd.h.
182 We do, however, guarantee that
183
184 Alen >= UMF_COLAMD_RECOMMENDED (nnz, n_row, n_col)
185
186 will be sufficient.
187
188 Int A [Alen] ; Input and output argument.
189
190 A is an integer array of size Alen. Alen must be at least as
191 large as the bare minimum value given above, but this is very
192 low, and can result in excessive run time. For best
193 performance, we recommend that Alen be greater than or equal to
194 UMF_COLAMD_RECOMMENDED (nnz, n_row, n_col), which adds
195 nnz/5 to the bare minimum value given above.
196
197 On input, the row indices of the entries in column c of the
198 matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices
199 in a given column c need not be in ascending order, and
200 duplicate row indices may be be present. However, colamd will
201 work a little faster if both of these conditions are met
202 (Colamd puts the matrix into this format, if it finds that the
203 the conditions are not met).
204
205 The matrix is 0-based. That is, rows are in the range 0 to
206 n_row-1, and columns are in the range 0 to n_col-1. Colamd
207 returns FALSE if any row index is out of range.
208
209 A holds the inverse permutation on output.
210
211 Int p [n_col+1] ; Both input and output argument.
212
213 p is an integer array of size n_col+1. On input, it holds the
214 "pointers" for the column form of the matrix A. Column c of
215 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
216 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
217 for all c in the range 0 to n_col-1. The value p [n_col] is
218 thus the total number of entries in the pattern of the matrix A.
219 Colamd returns FALSE if these conditions are not met.
220
221 On output, if colamd returns TRUE, the array p holds the column
222 permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
223 the first column index in the new ordering, and p [n_col-1] is
224 the last. That is, p [k] = j means that column j of A is the
225 kth pivot column, in AQ, where k is in the range 0 to n_col-1
226 (p [0] = j means that column j of A is the first column in AQ).
227
228 If colamd returns FALSE, then no permutation is returned, and
229 p is undefined on output.
230
231 double knobs [COLAMD_KNOBS] ; Input argument.
232
233 See colamd_set_defaults for a description.
234 The behavior is undefined if knobs contains NaN's.
235 (UMFPACK does not call umf_colamd with NaN-valued knobs).
236
237 Int stats [COLAMD_STATS] ; Output argument.
238
239 Statistics on the ordering, and error status.
240 See colamd.h for related definitions.
241 Colamd returns FALSE if stats is not present.
242
243 stats [0]: number of dense or empty rows ignored.
244
245 stats [1]: number of dense or empty columns ignored (and
246 ordered last in the output permutation p)
247 Note that a row can become "empty" if it
248 contains only "dense" and/or "empty" columns,
249 and similarly a column can become "empty" if it
250 only contains "dense" and/or "empty" rows.
251
252 stats [2]: number of garbage collections performed.
253 This can be excessively high if Alen is close
254 to the minimum required value.
255
256 stats [3]: status code. < 0 is an error code.
257 > 1 is a warning or notice.
258
259 0 OK. Each column of the input matrix contained
260 row indices in increasing order, with no
261 duplicates.
262
263 -11 Columns of input matrix jumbled
264 (unsorted columns or duplicate entries).
265
266 stats [4]: the bad column index
267 stats [5]: the bad row index
268
269 -1 A is a null pointer
270
271 -2 p is a null pointer
272
273 -3 n_row is negative
274
275 stats [4]: n_row
276
277 -4 n_col is negative
278
279 stats [4]: n_col
280
281 -5 number of nonzeros in matrix is negative
282
283 stats [4]: number of nonzeros, p [n_col]
284
285 -6 p [0] is nonzero
286
287 stats [4]: p [0]
288
289 -7 A is too small
290
291 stats [4]: required size
292 stats [5]: actual size (Alen)
293
294 -8 a column has a zero or negative number of
295 entries (changed for UMFPACK)
296
297 stats [4]: column with <= 0 entries
298 stats [5]: number of entries in col
299
300 -9 a row index is out of bounds
301
302 stats [4]: column with bad row index
303 stats [5]: bad row index
304 stats [6]: n_row, # of rows of matrx
305
306 -10 unused
307
308 -999 (unused; see symamd.c)
309
310 Future versions may return more statistics in the stats array.
311
312 Example:
313
314 See colamd_example.c for a complete example.
315
316 To order the columns of a 5-by-4 matrix with 11 nonzero entries in
317 the following nonzero pattern
318
319 x 0 x 0
320 x 0 x x
321 0 x x 0
322 0 0 x x
323 x x 0 0
324
325 with default knobs and no output statistics, do the following:
326
327 #include "colamd.h"
328 #define ALEN UMF_COLAMD_RECOMMENDED (11, 5, 4)
329 Int A [ALEN] = {1, 2, 5, 3, 5, 1, 2, 3, 4, 2, 4} ;
330 Int p [ ] = {0, 3, 5, 9, 11} ;
331 Int stats [COLAMD_STATS] ;
332 UMF_colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
333
334 The permutation is returned in the array p, and A is destroyed.
335
336
337 ----------------------------------------------------------------------------
338 symamd: does not appear in this version for UMFPACK
339 ----------------------------------------------------------------------------
340
341 ----------------------------------------------------------------------------
342 colamd_report: does not appear in this version for UMFPACK
343 ----------------------------------------------------------------------------
344
345 ----------------------------------------------------------------------------
346 symamd_report: does not appear in this version for UMFPACK
347 ----------------------------------------------------------------------------
348
349 */
350
351 /* ========================================================================== */
352 /* === Scaffolding code definitions ======================================== */
353 /* ========================================================================== */
354
355 /* UMFPACK debugging control moved to amd_internal.h */
356
357 /*
358 Our "scaffolding code" philosophy: In our opinion, well-written library
359 code should keep its "debugging" code, and just normally have it turned off
360 by the compiler so as not to interfere with performance. This serves
361 several purposes:
362
363 (1) assertions act as comments to the reader, telling you what the code
364 expects at that point. All assertions will always be true (unless
365 there really is a bug, of course).
366
367 (2) leaving in the scaffolding code assists anyone who would like to modify
368 the code, or understand the algorithm (by reading the debugging output,
369 one can get a glimpse into what the code is doing).
370
371 (3) (gasp!) for actually finding bugs. This code has been heavily tested
372 and "should" be fully functional and bug-free ... but you never know...
373
374 To enable debugging, comment out the "#define NDEBUG" above. For a MATLAB
375 mexFunction, you will also need to modify mexopts.sh to remove the -DNDEBUG
376 definition. The code will become outrageously slow when debugging is
377 enabled. To control the level of debugging output, set an environment
378 variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging,
379 you should see the following message on the standard output:
380
381 colamd: debug version, D = 1 (THIS WILL BE SLOW!)
382
383 or a similar message for symamd. If you don't, then debugging has not
384 been enabled.
385
386 */
387
388 /* ========================================================================== */
389 /* === Include files ======================================================== */
390 /* ========================================================================== */
391
392 /* ------------------ */
393 /* modified for UMFPACK: */
394 #include "umf_internal.h"
395 #include "umf_colamd.h"
396 #include "umf_apply_order.h"
397 #include "umf_fsize.h"
398 /* ------------------ */
399
400 /* ========================================================================== */
401 /* === Definitions ========================================================== */
402 /* ========================================================================== */
403
404 /* ------------------ */
405 /* UMFPACK: duplicate definitions moved to umf_internal.h */
406 /* ------------------ */
407
408 /* Row and column status */
409 #define ALIVE (0)
410 #define DEAD (-1)
411
412 /* Column status */
413 #define DEAD_PRINCIPAL (-1)
414 #define DEAD_NON_PRINCIPAL (-2)
415
416 /* Macros for row and column status update and checking. */
417 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
418 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
419 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
420 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
421 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
422 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
423 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
424 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
425 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
426
427 /* ------------------ */
428 /* UMFPACK: Colamd reporting mechanism moved to umf_internal.h */
429 /* ------------------ */
430
431 /* ========================================================================== */
432 /* === Prototypes of PRIVATE routines ======================================= */
433 /* ========================================================================== */
434
435 PRIVATE Int init_rows_cols
436 (
437 Int n_row,
438 Int n_col,
439 Colamd_Row Row [],
440 Colamd_Col Col [],
441 Int A [],
442 Int p []
443 /* Int stats [COLAMD_STATS] */
444 ) ;
445
446 PRIVATE void init_scoring
447 (
448 Int n_row,
449 Int n_col,
450 Colamd_Row Row [],
451 Colamd_Col Col [],
452 Int A [],
453 Int head [],
454 double knobs [COLAMD_KNOBS],
455 Int *p_n_row2,
456 Int *p_n_col2,
457 Int *p_max_deg
458 /* ------------------ */
459 /* added for UMFPACK */
460 , Int *p_ndense_row /* number of dense rows */
461 , Int *p_nempty_row /* number of original empty rows */
462 , Int *p_nnewlyempty_row /* number of newly empty rows */
463 , Int *p_ndense_col /* number of dense cols (excl "empty" cols) */
464 , Int *p_nempty_col /* number of original empty cols */
465 , Int *p_nnewlyempty_col /* number of newly empty cols */
466 ) ;
467
468 PRIVATE Int find_ordering
469 (
470 Int n_row,
471 Int n_col,
472 Int Alen,
473 Colamd_Row Row [],
474 Colamd_Col Col [],
475 Int A [],
476 Int head [],
477 Int n_col2,
478 Int max_deg,
479 Int pfree
480 /* ------------------ */
481 /* added for UMFPACK: */
482 , Int Front_npivcol [ ]
483 , Int Front_nrows [ ]
484 , Int Front_ncols [ ]
485 , Int Front_parent [ ]
486 , Int Front_cols [ ]
487 , Int *p_nfr
488 , Int aggressive
489 , Int InFront [ ]
490 /* ------------------ */
491 ) ;
492
493 /* ------------------ */
494 /* order_children deleted for UMFPACK: */
495 /* ------------------ */
496
497 PRIVATE void detect_super_cols
498 (
499
500 #ifndef NDEBUG
501 Int n_col,
502 Colamd_Row Row [],
503 #endif /* NDEBUG */
504
505 Colamd_Col Col [],
506 Int A [],
507 Int head [],
508 Int row_start,
509 Int row_length
510 ) ;
511
512 PRIVATE Int garbage_collection
513 (
514 Int n_row,
515 Int n_col,
516 Colamd_Row Row [],
517 Colamd_Col Col [],
518 Int A [],
519 Int *pfree
520 ) ;
521
522 PRIVATE Int clear_mark
523 (
524 Int n_row,
525 Colamd_Row Row []
526 ) ;
527
528 /* ------------------ */
529 /* print_report deleted for UMFPACK */
530 /* ------------------ */
531
532 /* ========================================================================== */
533 /* === Debugging prototypes and definitions ================================= */
534 /* ========================================================================== */
535
536 #ifndef NDEBUG
537
538 /* ------------------ */
539 /* debugging macros moved for UMFPACK */
540 /* ------------------ */
541
542 PRIVATE void debug_deg_lists
543 (
544 Int n_row,
545 Int n_col,
546 Colamd_Row Row [],
547 Colamd_Col Col [],
548 Int head [],
549 Int min_score,
550 Int should,
551 Int max_deg
552 ) ;
553
554 PRIVATE void debug_mark
555 (
556 Int n_row,
557 Colamd_Row Row [],
558 Int tag_mark,
559 Int max_mark
560 ) ;
561
562 PRIVATE void debug_matrix
563 (
564 Int n_row,
565 Int n_col,
566 Colamd_Row Row [],
567 Colamd_Col Col [],
568 Int A []
569 ) ;
570
571 PRIVATE void debug_structures
572 (
573 Int n_row,
574 Int n_col,
575 Colamd_Row Row [],
576 Colamd_Col Col [],
577 Int A [],
578 Int n_col2
579 ) ;
580
581 /* ------------------ */
582 /* dump_super added for UMFPACK: */
583 PRIVATE void dump_super
584 (
585 Int super_c,
586 Colamd_Col Col [],
587 Int n_col
588 ) ;
589 /* ------------------ */
590
591 #endif /* NDEBUG */
592
593 /* ========================================================================== */
594
595
596
597 /* ========================================================================== */
598 /* === USER-CALLABLE ROUTINES: ============================================== */
599 /* ========================================================================== */
600
601
602 /* ========================================================================== */
603 /* === colamd_set_defaults ================================================== */
604 /* ========================================================================== */
605
606 /*
607 The colamd_set_defaults routine sets the default values of the user-
608 controllable parameters for colamd:
609
610 knobs [0] rows with knobs[0]*n_col entries or more are removed
611 prior to ordering in colamd. Rows and columns with
612 knobs[0]*n_col entries or more are removed prior to
613 ordering in symamd and placed last in the output
614 ordering.
615
616 knobs [1] columns with knobs[1]*n_row entries or more are removed
617 prior to ordering in colamd, and placed last in the
618 column permutation. Symamd ignores this knob.
619
620 knobs [2] if nonzero, then perform aggressive absorption.
621
622 knobs [3..19] unused, but future versions might use this
623 */
624
UMF_colamd_set_defaults(double knobs[COLAMD_KNOBS])625 GLOBAL void UMF_colamd_set_defaults
626 (
627 /* === Parameters ======================================================= */
628
629 double knobs [COLAMD_KNOBS] /* knob array */
630 )
631 {
632 /* === Local variables ================================================== */
633
634 Int i ;
635
636 #if 0
637 if (!knobs)
638 {
639 return ; /* UMFPACK always passes knobs array */
640 }
641 #endif
642 for (i = 0 ; i < COLAMD_KNOBS ; i++)
643 {
644 knobs [i] = 0 ;
645 }
646 knobs [COLAMD_DENSE_ROW] = 0.2 ; /* default changed for UMFPACK */
647 knobs [COLAMD_DENSE_COL] = 0.2 ; /* default changed for UMFPACK */
648 knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default is to do aggressive
649 * absorption */
650 }
651
652
653 /* ========================================================================== */
654 /* === symamd removed for UMFPACK =========================================== */
655 /* ========================================================================== */
656
657
658
659 /* ========================================================================== */
660 /* === colamd =============================================================== */
661 /* ========================================================================== */
662
663 /*
664 The colamd routine computes a column ordering Q of a sparse matrix
665 A such that the LU factorization P(AQ) = LU remains sparse, where P is
666 selected via partial pivoting. The routine can also be viewed as
667 providing a permutation Q such that the Cholesky factorization
668 (AQ)'(AQ) = LL' remains sparse.
669 */
670
671 /* For UMFPACK: colamd always returns TRUE */
672
UMF_colamd(Int n_row,Int n_col,Int Alen,Int A[],Int p[],double knobs[COLAMD_KNOBS],Int stats[COLAMD_STATS],Int Front_npivcol[],Int Front_nrows[],Int Front_ncols[],Int Front_parent[],Int Front_cols[],Int * p_nfr,Int InFront[])673 GLOBAL Int UMF_colamd /* returns TRUE if successful, FALSE otherwise*/
674 (
675 /* === Parameters ======================================================= */
676
677 Int n_row, /* number of rows in A */
678 Int n_col, /* number of columns in A */
679 Int Alen, /* length of A */
680 Int A [], /* row indices of A */
681 Int p [], /* pointers to columns in A */
682 double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
683 Int stats [COLAMD_STATS] /* output statistics and error codes */
684
685 /* ------------------ */
686 /* added for UMFPACK: each Front_ array is of size n_col+1 */
687 , Int Front_npivcol [ ] /* # pivot cols in each front */
688 , Int Front_nrows [ ] /* # of rows in each front (incl. pivot rows) */
689 , Int Front_ncols [ ] /* # of cols in each front (incl. pivot cols) */
690 , Int Front_parent [ ] /* parent of each front */
691 , Int Front_cols [ ] /* link list of pivot columns for each front */
692 , Int *p_nfr /* total number of frontal matrices */
693 , Int InFront [ ] /* InFront [row] = f if the original row was
694 * absorbed into front f. EMPTY if the row was
695 * empty, dense, or not absorbed. This array
696 * has size n_row+1 */
697 /* ------------------ */
698 )
699 {
700 /* === Local variables ================================================== */
701
702 Int row ; /* row index */
703 Int i ; /* loop index */
704 Int nnz ; /* nonzeros in A */
705 Int Row_size ; /* size of Row [], in integers */
706 Int Col_size ; /* size of Col [], in integers */
707 #if 0
708 Int need ; /* minimum required length of A */
709 #endif
710 Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */
711 Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */
712 Int n_col2 ; /* number of non-dense, non-empty columns */
713 Int n_row2 ; /* number of non-dense, non-empty rows */
714 Int ngarbage ; /* number of garbage collections performed */
715 Int max_deg ; /* maximum row degree */
716 Int aggressive ; /* TRUE if doing aggressive absorption */
717 #if 0
718 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
719 #endif
720
721 /* ------------------ */
722 /* debugging initializations moved for UMFPACK */
723 /* ------------------ */
724
725 /* ------------------ */
726 /* added for UMFPACK: */
727 Int ndense_row, nempty_row, parent, ndense_col,
728 nempty_col, k, col, nfr, *Front_child, *Front_sibling, *Front_stack,
729 *Front_order, *Front_size ;
730 Int nnewlyempty_col, nnewlyempty_row ;
731 /* ------------------ */
732
733 /* === Check the input arguments ======================================== */
734
735 #if 0
736 if (!stats)
737 {
738 DEBUG0 (("colamd: stats not present\n")) ;
739 return (FALSE) ; /* UMFPACK: always passes stats [ ] */
740 }
741 #endif
742
743 ASSERT (stats != (Int *) NULL) ;
744
745 for (i = 0 ; i < COLAMD_STATS ; i++)
746 {
747 stats [i] = 0 ;
748 }
749 stats [COLAMD_STATUS] = COLAMD_OK ;
750 stats [COLAMD_INFO1] = -1 ;
751 stats [COLAMD_INFO2] = -1 ;
752
753 #if 0
754 if (!A) /* A is not present */
755 {
756 /* UMFPACK: always passes A [ ] */
757 DEBUG0 (("colamd: A not present\n")) ;
758 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
759 return (FALSE) ;
760 }
761
762 if (!p) /* p is not present */
763 {
764 /* UMFPACK: always passes p [ ] */
765 DEBUG0 (("colamd: p not present\n")) ;
766 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
767 return (FALSE) ;
768 }
769
770 if (n_row < 0) /* n_row must be >= 0 */
771 {
772 /* UMFPACK: does not call UMF_colamd if n <= 0 */
773 DEBUG0 (("colamd: nrow negative "ID"\n", n_row)) ;
774 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
775 stats [COLAMD_INFO1] = n_row ;
776 return (FALSE) ;
777 }
778
779 if (n_col < 0) /* n_col must be >= 0 */
780 {
781 /* UMFPACK: does not call UMF_colamd if n <= 0 */
782 DEBUG0 (("colamd: ncol negative "ID"\n", n_col)) ;
783 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
784 stats [COLAMD_INFO1] = n_col ;
785 return (FALSE) ;
786 }
787 #endif
788
789 ASSERT (A != (Int *) NULL) ;
790 ASSERT (p != (Int *) NULL) ;
791 ASSERT (n_row >= 0) ;
792 ASSERT (n_col >= 0) ;
793
794 nnz = p [n_col] ;
795
796 #if 0
797 if (nnz < 0) /* nnz must be >= 0 */
798 {
799 /* UMFPACK: does not call UMF_colamd if nnz < 0 */
800 DEBUG0 (("colamd: number of entries negative "ID"\n", nnz)) ;
801 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
802 stats [COLAMD_INFO1] = nnz ;
803 return (FALSE) ;
804 }
805
806 if (p [0] != 0) /* p [0] must be exactly zero */
807 {
808 DEBUG0 (("colamd: p[0] not zero "ID"\n", p [0])) ;
809 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
810 stats [COLAMD_INFO1] = p [0] ;
811 return (FALSE) ;
812 }
813 #endif
814
815 ASSERT (nnz >= 0) ;
816 ASSERT (p [0] == 0) ;
817
818 /* === If no knobs, set default knobs =================================== */
819
820 #if 0
821 if (!knobs)
822 {
823 /* UMFPACK: always passes the knobs */
824 UMF_colamd_set_defaults (default_knobs) ;
825 knobs = default_knobs ;
826 }
827 #endif
828
829 ASSERT (knobs != (double *) NULL) ;
830
831 /* --------------------- */
832 /* added for UMFPACK v4.1: */
833 aggressive = (knobs [COLAMD_AGGRESSIVE] != 0) ;
834 /* --------------------- */
835
836 /* === Allocate the Row and Col arrays from array A ===================== */
837
838 Col_size = UMF_COLAMD_C (n_col) ;
839 Row_size = UMF_COLAMD_R (n_row) ;
840
841 #if 0
842 need = MAX (2*nnz, 4*n_col) + n_col + Col_size + Row_size ;
843 if (need > Alen)
844 {
845 /* UMFPACK: always passes enough space */
846 /* not enough space in array A to perform the ordering */
847 DEBUG0 (("colamd: Need Alen >= "ID", given only Alen = "ID"\n",
848 need, Alen)) ;
849 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
850 stats [COLAMD_INFO1] = need ;
851 stats [COLAMD_INFO2] = Alen ;
852 return (FALSE) ;
853 }
854 #endif
855
856 Alen -= Col_size + Row_size ;
857 Col = (Colamd_Col *) &A [Alen] ;
858 Row = (Colamd_Row *) &A [Alen + Col_size] ;
859
860 /* Size of A is now Alen >= MAX (2*nnz, 4*n_col) + n_col. The ordering
861 * requires Alen >= 2*nnz + n_col, and the postorder requires
862 * Alen >= 5*n_col. */
863
864 /* === Construct the row and column data structures ===================== */
865
866 i = init_rows_cols (n_row, n_col, Row, Col, A, p) ;
867
868 #if 0
869 if (!i)
870 {
871 /* input matrix is invalid */
872 DEBUG0 (("colamd: Matrix invalid\n")) ;
873 return (FALSE) ;
874 }
875 #endif
876
877 ASSERT (i) ;
878
879 /* === UMFPACK: Initialize front info =================================== */
880
881 for (col = 0 ; col < n_col ; col++)
882 {
883 Front_npivcol [col] = 0 ;
884 Front_nrows [col] = 0 ;
885 Front_ncols [col] = 0 ;
886 Front_parent [col] = EMPTY ;
887 Front_cols [col] = EMPTY ;
888 }
889
890 /* === Initialize scores, kill dense rows/columns ======================= */
891
892 init_scoring (n_row, n_col, Row, Col, A, p, knobs,
893 &n_row2, &n_col2, &max_deg
894 /* ------------------ */
895 /* added for UMFPACK: */
896 , &ndense_row, &nempty_row, &nnewlyempty_row
897 , &ndense_col, &nempty_col, &nnewlyempty_col
898 /* ------------------ */
899 ) ;
900 ASSERT (n_row2 == n_row - nempty_row - nnewlyempty_row - ndense_row) ;
901 ASSERT (n_col2 == n_col - nempty_col - nnewlyempty_col - ndense_col) ;
902
903 /* === Order the supercolumns =========================================== */
904
905 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
906 n_col2, max_deg, 2*nnz
907 /* ------------------ */
908 /* added for UMFPACK: */
909 , Front_npivcol, Front_nrows, Front_ncols, Front_parent, Front_cols
910 , &nfr, aggressive, InFront
911 /* ------------------ */
912 ) ;
913
914 /* ------------------ */
915 /* changed for UMFPACK: */
916
917 /* A is no longer needed, so use A [0..5*nfr-1] as workspace [ [ */
918 /* This step requires Alen >= 5*n_col */
919 Front_child = A ;
920 Front_sibling = Front_child + nfr ;
921 Front_stack = Front_sibling + nfr ;
922 Front_order = Front_stack + nfr ;
923 Front_size = Front_order + nfr ;
924
925 UMF_fsize (nfr, Front_size, Front_nrows, Front_ncols,
926 Front_parent, Front_npivcol) ;
927
928 AMD_postorder (nfr, Front_parent, Front_npivcol, Front_size,
929 Front_order, Front_child, Front_sibling, Front_stack) ;
930
931 /* Front_size, Front_stack, Front_child, Front_sibling no longer needed ] */
932
933 /* use A [0..nfr-1] as workspace */
934 UMF_apply_order (Front_npivcol, Front_order, A, nfr, nfr) ;
935 UMF_apply_order (Front_nrows, Front_order, A, nfr, nfr) ;
936 UMF_apply_order (Front_ncols, Front_order, A, nfr, nfr) ;
937 UMF_apply_order (Front_parent, Front_order, A, nfr, nfr) ;
938 UMF_apply_order (Front_cols, Front_order, A, nfr, nfr) ;
939
940 /* fix the parent to refer to the new numbering */
941 for (i = 0 ; i < nfr ; i++)
942 {
943 parent = Front_parent [i] ;
944 if (parent != EMPTY)
945 {
946 Front_parent [i] = Front_order [parent] ;
947 }
948 }
949
950 /* fix InFront to refer to the new numbering */
951 for (row = 0 ; row < n_row ; row++)
952 {
953 i = InFront [row] ;
954 ASSERT (i >= EMPTY && i < nfr) ;
955 if (i != EMPTY)
956 {
957 InFront [row] = Front_order [i] ;
958 }
959 }
960
961 /* Front_order longer needed ] */
962
963 /* === Order the columns in the fronts ================================== */
964
965 /* use A [0..n_col-1] as inverse permutation */
966 for (i = 0 ; i < n_col ; i++)
967 {
968 A [i] = EMPTY ;
969 }
970 k = 0 ;
971 for (i = 0 ; i < nfr ; i++)
972 {
973 ASSERT (Front_npivcol [i] > 0) ;
974 for (col = Front_cols [i] ; col != EMPTY ; col = Col [col].nextcol)
975 {
976 ASSERT (col >= 0 && col < n_col) ;
977 DEBUG1 (("Colamd output ordering: k "ID" col "ID"\n", k, col)) ;
978 p [k] = col ;
979 ASSERT (A [col] == EMPTY) ;
980 A [col] = k ;
981 k++ ;
982 }
983 }
984
985 /* === Order the "dense" and null columns =============================== */
986
987 ASSERT (k == n_col2) ;
988 if (n_col2 < n_col)
989 {
990 for (col = 0 ; col < n_col ; col++)
991 {
992 if (A [col] == EMPTY)
993 {
994 k = Col [col].shared2.order ;
995 ASSERT (k >= n_col2 && k < n_col) ;
996 DEBUG1 (("Colamd output ordering: k "ID" col "ID
997 " (dense or null col)\n", k, col)) ;
998 p [k] = col ;
999 A [col] = k ;
1000 }
1001 }
1002 }
1003
1004 /* ------------------ */
1005
1006 /* === Return statistics in stats ======================================= */
1007
1008 /* ------------------ */
1009 /* modified for UMFPACK */
1010 stats [COLAMD_DENSE_ROW] = ndense_row ;
1011 stats [COLAMD_EMPTY_ROW] = nempty_row ;
1012 stats [COLAMD_NEWLY_EMPTY_ROW] = nnewlyempty_row ;
1013 stats [COLAMD_DENSE_COL] = ndense_col ;
1014 stats [COLAMD_EMPTY_COL] = nempty_col ;
1015 stats [COLAMD_NEWLY_EMPTY_COL] = nnewlyempty_col ;
1016 ASSERT (ndense_col + nempty_col + nnewlyempty_col == n_col - n_col2) ;
1017 /* ------------------ */
1018 stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
1019 *p_nfr = nfr ;
1020 DEBUG1 (("colamd: done.\n")) ;
1021 return (TRUE) ;
1022 }
1023
1024
1025
1026
1027 /* ========================================================================== */
1028 /* === colamd_report removed for UMFPACK ==================================== */
1029 /* ========================================================================== */
1030
1031 /* ========================================================================== */
1032 /* === symamd_report removed for UMFPACK ==================================== */
1033 /* ========================================================================== */
1034
1035
1036
1037 /* ========================================================================== */
1038 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
1039 /* ========================================================================== */
1040
1041 /* There are no user-callable routines beyond this point in the file */
1042
1043
1044 /* ========================================================================== */
1045 /* === init_rows_cols ======================================================= */
1046 /* ========================================================================== */
1047
1048 /*
1049 Takes the column form of the matrix in A and creates the row form of the
1050 matrix. Also, row and column attributes are stored in the Col and Row
1051 structs. If the columns are un-sorted or contain duplicate row indices,
1052 this routine will also sort and remove duplicate row indices from the
1053 column form of the matrix. Returns FALSE if the matrix is invalid,
1054 TRUE otherwise. Not user-callable.
1055 */
1056
1057 /* For UMFPACK, this always returns TRUE */
1058
init_rows_cols(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int p[])1059 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
1060 (
1061 /* === Parameters ======================================================= */
1062
1063 Int n_row, /* number of rows of A */
1064 Int n_col, /* number of columns of A */
1065 Colamd_Row Row [], /* of size n_row+1 */
1066 Colamd_Col Col [], /* of size n_col+1 */
1067 Int A [], /* row indices of A, of size Alen */
1068 Int p [] /* pointers to columns in A, of size n_col+1 */
1069 /*
1070 Int stats [COLAMD_STATS] colamd statistics, removed for UMFPACK
1071 */
1072 )
1073 {
1074 /* === Local variables ================================================== */
1075
1076 Int col ; /* a column index */
1077 Int row ; /* a row index */
1078 Int *cp ; /* a column pointer */
1079 Int *cp_end ; /* a pointer to the end of a column */
1080
1081 /* === Initialize columns, and check column pointers ==================== */
1082
1083 for (col = 0 ; col < n_col ; col++)
1084 {
1085 Col [col].start = p [col] ;
1086 Col [col].length = p [col+1] - p [col] ;
1087
1088 #if 0
1089 if (Col [col].length < 0)
1090 {
1091 /* column pointers must be non-decreasing */
1092 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1093 stats [COLAMD_INFO1] = col ;
1094 stats [COLAMD_INFO2] = Col [col].length ;
1095 DEBUG0 (("colamd: col "ID" length "ID" <= 0\n",
1096 col, Col [col].length));
1097 return (FALSE) ;
1098 }
1099 #endif
1100
1101 ASSERT (Col [col].length >= 0) ;
1102
1103 /* added for UMFPACK v4.1 */
1104 ASSERT (Col [col].length > 0) ;
1105
1106 Col [col].shared1.thickness = 1 ;
1107 Col [col].shared2.score = 0 ;
1108 Col [col].shared3.prev = EMPTY ;
1109 Col [col].shared4.degree_next = EMPTY ;
1110
1111 /* ------------------ */
1112 /* added for UMFPACK: */
1113 Col [col].nextcol = EMPTY ;
1114 Col [col].lastcol = col ;
1115 /* ------------------ */
1116 }
1117
1118 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
1119
1120 /* === Scan columns, compute row degrees, and check row indices ========= */
1121
1122 /* ------------------ */
1123 /* stats [COLAMD_INFO3] = 0 ; */
1124 /* number of duplicate or unsorted row indices - not computed in UMFPACK */
1125 /* ------------------ */
1126
1127 for (row = 0 ; row < n_row ; row++)
1128 {
1129 Row [row].length = 0 ;
1130 /* ------------------ */
1131 /* removed for UMFPACK */
1132 /* Row [row].shared2.mark = -1 ; */
1133 /* ------------------ */
1134 /* ------------------ */
1135 /* added for UMFPACK: */
1136 Row [row].thickness = 1 ;
1137 Row [row].front = EMPTY ;
1138 /* ------------------ */
1139 }
1140
1141 for (col = 0 ; col < n_col ; col++)
1142 {
1143 #ifndef NDEBUG
1144 Int last_row = -1 ;
1145 #endif
1146
1147 cp = &A [p [col]] ;
1148 cp_end = &A [p [col+1]] ;
1149
1150 while (cp < cp_end)
1151 {
1152 row = *cp++ ;
1153
1154 #if 0
1155 /* make sure row indices within range */
1156 if (row < 0 || row >= n_row)
1157 {
1158 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1159 stats [COLAMD_INFO1] = col ;
1160 stats [COLAMD_INFO2] = row ;
1161 /* ------------------ */
1162 /* not needed in UMFPACK: */
1163 /* stats [COLAMD_INFO3] = n_row ; */
1164 /* ------------------ */
1165 DEBUG0 (("colamd: row "ID" col "ID" out of bounds\n", row,col));
1166 return (FALSE) ;
1167 }
1168 #endif
1169
1170 ASSERT (row >= 0 && row < n_row) ;
1171
1172 #if 0
1173 /* ------------------ */
1174 /* changed for UMFPACK */
1175 if (row <= last_row)
1176 {
1177 /* row index are unsorted or repeated (or both), thus col */
1178 /* is jumbled. This is an error condition for UMFPACK */
1179 stats [COLAMD_STATUS] = COLAMD_ERROR_jumbled_matrix ;
1180 stats [COLAMD_INFO1] = col ;
1181 stats [COLAMD_INFO2] = row ;
1182 DEBUG1 (("colamd: row "ID" col "ID" unsorted/duplicate\n",
1183 row, col)) ;
1184 return (FALSE) ;
1185 }
1186 /* ------------------ */
1187 #endif
1188
1189 ASSERT (row > last_row) ;
1190
1191 /* ------------------ */
1192 /* changed for UMFPACK - jumbled columns not tolerated */
1193 Row [row].length++ ;
1194 /* ------------------ */
1195
1196 #ifndef NDEBUG
1197 last_row = row ;
1198 #endif
1199 }
1200 }
1201
1202 /* === Compute row pointers ============================================= */
1203
1204 /* row form of the matrix starts directly after the column */
1205 /* form of matrix in A */
1206 Row [0].start = p [n_col] ;
1207 Row [0].shared1.p = Row [0].start ;
1208 /* ------------------ */
1209 /* removed for UMFPACK */
1210 /* Row [0].shared2.mark = -1 ; */
1211 /* ------------------ */
1212 for (row = 1 ; row < n_row ; row++)
1213 {
1214 Row [row].start = Row [row-1].start + Row [row-1].length ;
1215 Row [row].shared1.p = Row [row].start ;
1216 /* ------------------ */
1217 /* removed for UMFPACK */
1218 /* Row [row].shared2.mark = -1 ; */
1219 /* ------------------ */
1220 }
1221
1222 /* === Create row form ================================================== */
1223
1224 /* ------------------ */
1225 /* jumbled matrix case removed for UMFPACK */
1226 /* ------------------ */
1227
1228 for (col = 0 ; col < n_col ; col++)
1229 {
1230 cp = &A [p [col]] ;
1231 cp_end = &A [p [col+1]] ;
1232 while (cp < cp_end)
1233 {
1234 A [(Row [*cp++].shared1.p)++] = col ;
1235 }
1236 }
1237
1238 /* === Clear the row marks and set row degrees ========================== */
1239
1240 for (row = 0 ; row < n_row ; row++)
1241 {
1242 Row [row].shared2.mark = 0 ;
1243 Row [row].shared1.degree = Row [row].length ;
1244 }
1245
1246 /* ------------------ */
1247 /* recreate columns for jumbled matrix case removed for UMFPACK */
1248 /* ------------------ */
1249
1250 return (TRUE) ;
1251 }
1252
1253
1254 /* ========================================================================== */
1255 /* === init_scoring ========================================================= */
1256 /* ========================================================================== */
1257
1258 /*
1259 Kills dense or empty columns and rows, calculates an initial score for
1260 each column, and places all columns in the degree lists. Not user-callable.
1261 */
1262
init_scoring(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int head[],double knobs[COLAMD_KNOBS],Int * p_n_row2,Int * p_n_col2,Int * p_max_deg,Int * p_ndense_row,Int * p_nempty_row,Int * p_nnewlyempty_row,Int * p_ndense_col,Int * p_nempty_col,Int * p_nnewlyempty_col)1263 PRIVATE void init_scoring
1264 (
1265 /* === Parameters ======================================================= */
1266
1267 Int n_row, /* number of rows of A */
1268 Int n_col, /* number of columns of A */
1269 Colamd_Row Row [], /* of size n_row+1 */
1270 Colamd_Col Col [], /* of size n_col+1 */
1271 Int A [], /* column form and row form of A */
1272 Int head [], /* of size n_col+1 */
1273 double knobs [COLAMD_KNOBS],/* parameters */
1274 Int *p_n_row2, /* number of non-dense, non-empty rows */
1275 Int *p_n_col2, /* number of non-dense, non-empty columns */
1276 Int *p_max_deg /* maximum row degree */
1277 /* ------------------ */
1278 /* added for UMFPACK */
1279 , Int *p_ndense_row /* number of dense rows */
1280 , Int *p_nempty_row /* number of original empty rows */
1281 , Int *p_nnewlyempty_row /* number of newly empty rows */
1282 , Int *p_ndense_col /* number of dense cols (excl "empty" cols) */
1283 , Int *p_nempty_col /* number of original empty cols */
1284 , Int *p_nnewlyempty_col /* number of newly empty cols */
1285 /* ------------------ */
1286 )
1287 {
1288 /* === Local variables ================================================== */
1289
1290 Int c ; /* a column index */
1291 Int r, row ; /* a row index */
1292 Int *cp ; /* a column pointer */
1293 Int deg ; /* degree of a row or column */
1294 Int *cp_end ; /* a pointer to the end of a column */
1295 Int *new_cp ; /* new column pointer */
1296 Int col_length ; /* length of pruned column */
1297 Int score ; /* current column score */
1298 Int n_col2 ; /* number of non-dense, non-empty columns */
1299 Int n_row2 ; /* number of non-dense, non-empty rows */
1300 Int dense_row_count ; /* remove rows with more entries than this */
1301 Int dense_col_count ; /* remove cols with more entries than this */
1302 Int min_score ; /* smallest column score */
1303 Int max_deg ; /* maximum row degree */
1304 Int next_col ; /* Used to add to degree list.*/
1305
1306 /* ------------------ */
1307 /* added for UMFPACK */
1308 Int ndense_row ; /* number of dense rows */
1309 Int nempty_row ; /* number of empty rows */
1310 Int nnewlyempty_row ; /* number of newly empty rows */
1311 Int ndense_col ; /* number of dense cols (excl "empty" cols) */
1312 Int nempty_col ; /* number of original empty cols */
1313 Int nnewlyempty_col ; /* number of newly empty cols */
1314 Int ne ;
1315 /* ------------------ */
1316
1317 #ifndef NDEBUG
1318 Int debug_count ; /* debug only. */
1319 #endif /* NDEBUG */
1320
1321 /* === Extract knobs ==================================================== */
1322
1323 /* --------------------- */
1324 /* old dense row/column knobs:
1325 dense_row_count = MAX (0, MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ;
1326 dense_col_count = MAX (0, MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ;
1327 */
1328 /* new, for UMFPACK: */
1329 /* Note: if knobs contains a NaN, this is undefined: */
1330 dense_row_count =
1331 UMFPACK_DENSE_DEGREE_THRESHOLD (knobs [COLAMD_DENSE_ROW], n_col) ;
1332 dense_col_count =
1333 UMFPACK_DENSE_DEGREE_THRESHOLD (knobs [COLAMD_DENSE_COL], n_row) ;
1334 /* Make sure dense_*_count is between 0 and n: */
1335 dense_row_count = MAX (0, MIN (dense_row_count, n_col)) ;
1336 dense_col_count = MAX (0, MIN (dense_col_count, n_row)) ;
1337 /* --------------------- */
1338
1339 DEBUG1 (("colamd: densecount: "ID" "ID"\n",
1340 dense_row_count, dense_col_count)) ;
1341 max_deg = 0 ;
1342 n_col2 = n_col ;
1343 n_row2 = n_row ;
1344
1345 /* --------------------- */
1346 /* added for UMFPACK */
1347 ndense_col = 0 ;
1348 nempty_col = 0 ;
1349 nnewlyempty_col = 0 ;
1350 ndense_row = 0 ;
1351 nempty_row = 0 ;
1352 nnewlyempty_row = 0 ;
1353 /* --------------------- */
1354
1355 /* === Kill empty columns =============================================== */
1356
1357 /* removed for UMFPACK v4.1. prune_singletons has already removed empty
1358 * columns and empty rows */
1359
1360 #if 0
1361 /* Put the empty columns at the end in their natural order, so that LU */
1362 /* factorization can proceed as far as possible. */
1363 for (c = n_col-1 ; c >= 0 ; c--)
1364 {
1365 deg = Col [c].length ;
1366 if (deg == 0)
1367 {
1368 /* this is a empty column, kill and order it last */
1369 Col [c].shared2.order = --n_col2 ;
1370 KILL_PRINCIPAL_COL (c) ;
1371 /* --------------------- */
1372 /* added for UMFPACK */
1373 nempty_col++ ;
1374 /* --------------------- */
1375 }
1376 }
1377 DEBUG1 (("colamd: null columns killed: "ID"\n", n_col - n_col2)) ;
1378 #endif
1379
1380 #ifndef NDEBUG
1381 for (c = 0 ; c < n_col ; c++)
1382 {
1383 ASSERT (Col [c].length > 0) ;
1384 }
1385 #endif
1386
1387 /* === Count null rows ================================================== */
1388
1389 #if 0
1390 for (r = 0 ; r < n_row ; r++)
1391 {
1392 deg = Row [r].shared1.degree ;
1393 if (deg == 0)
1394 {
1395 /* this is an original empty row */
1396 nempty_row++ ;
1397 }
1398 }
1399 #endif
1400
1401 #ifndef NDEBUG
1402 for (r = 0 ; r < n_row ; r++)
1403 {
1404 ASSERT (Row [r].shared1.degree > 0) ;
1405 ASSERT (Row [r].length > 0) ;
1406 }
1407 #endif
1408
1409 /* === Kill dense columns =============================================== */
1410
1411 /* Put the dense columns at the end, in their natural order */
1412 for (c = n_col-1 ; c >= 0 ; c--)
1413 {
1414
1415 /* ----------------------------------------------------------------- */
1416 #if 0
1417 /* removed for UMFPACK v4.1: no empty columns */
1418 /* skip any dead columns */
1419 if (COL_IS_DEAD (c))
1420 {
1421 continue ;
1422 }
1423 #endif
1424 ASSERT (COL_IS_ALIVE (c)) ;
1425 ASSERT (Col [c].length > 0) ;
1426 /* ----------------------------------------------------------------- */
1427
1428 deg = Col [c].length ;
1429 if (deg > dense_col_count)
1430 {
1431 /* this is a dense column, kill and order it last */
1432 Col [c].shared2.order = --n_col2 ;
1433 /* --------------------- */
1434 /* added for UMFPACK */
1435 ndense_col++ ;
1436 /* --------------------- */
1437 /* decrement the row degrees */
1438 cp = &A [Col [c].start] ;
1439 cp_end = cp + Col [c].length ;
1440 while (cp < cp_end)
1441 {
1442 Row [*cp++].shared1.degree-- ;
1443 }
1444 KILL_PRINCIPAL_COL (c) ;
1445 }
1446 }
1447 DEBUG1 (("colamd: Dense and null columns killed: "ID"\n", n_col - n_col2)) ;
1448
1449 /* === Kill dense and empty rows ======================================== */
1450
1451 /* Note that there can now be empty rows, since dense columns have
1452 * been deleted. These are "newly" empty rows. */
1453
1454 ne = 0 ;
1455 for (r = 0 ; r < n_row ; r++)
1456 {
1457 deg = Row [r].shared1.degree ;
1458 ASSERT (deg >= 0 && deg <= n_col) ;
1459 /* --------------------- */
1460 /* added for UMFPACK */
1461 if (deg > dense_row_count)
1462 {
1463 /* There is at least one dense row. Continue ordering, but */
1464 /* symbolic factorization will be redone after UMF_colamd is done.*/
1465 ndense_row++ ;
1466 }
1467 if (deg == 0)
1468 {
1469 /* this is a newly empty row, or original empty row */
1470 ne++ ;
1471 }
1472 /* --------------------- */
1473 if (deg > dense_row_count || deg == 0)
1474 {
1475 /* kill a dense or empty row */
1476 KILL_ROW (r) ;
1477 /* --------------------- */
1478 /* added for UMFPACK */
1479 Row [r].thickness = 0 ;
1480 /* --------------------- */
1481 --n_row2 ;
1482 }
1483 else
1484 {
1485 /* keep track of max degree of remaining rows */
1486 max_deg = MAX (max_deg, deg) ;
1487 }
1488 }
1489 nnewlyempty_row = ne - nempty_row ;
1490 DEBUG1 (("colamd: Dense rows killed: "ID"\n", ndense_row)) ;
1491 DEBUG1 (("colamd: Dense and null rows killed: "ID"\n", n_row - n_row2)) ;
1492
1493 /* === Compute initial column scores ==================================== */
1494
1495 /* At this point the row degrees are accurate. They reflect the number */
1496 /* of "live" (non-dense) columns in each row. No empty rows exist. */
1497 /* Some "live" columns may contain only dead rows, however. These are */
1498 /* pruned in the code below. */
1499
1500 /* now find the initial matlab score for each column */
1501 for (c = n_col-1 ; c >= 0 ; c--)
1502 {
1503 /* skip dead column */
1504 if (COL_IS_DEAD (c))
1505 {
1506 continue ;
1507 }
1508 score = 0 ;
1509 cp = &A [Col [c].start] ;
1510 new_cp = cp ;
1511 cp_end = cp + Col [c].length ;
1512 while (cp < cp_end)
1513 {
1514 /* get a row */
1515 row = *cp++ ;
1516 /* skip if dead */
1517 if (ROW_IS_DEAD (row))
1518 {
1519 continue ;
1520 }
1521 /* compact the column */
1522 *new_cp++ = row ;
1523 /* add row's external degree */
1524 score += Row [row].shared1.degree - 1 ;
1525 /* guard against integer overflow */
1526 score = MIN (score, n_col) ;
1527 }
1528 /* determine pruned column length */
1529 col_length = (Int) (new_cp - &A [Col [c].start]) ;
1530 if (col_length == 0)
1531 {
1532 /* a newly-made null column (all rows in this col are "dense" */
1533 /* and have already been killed) */
1534 DEBUG2 (("Newly null killed: "ID"\n", c)) ;
1535 Col [c].shared2.order = --n_col2 ;
1536 KILL_PRINCIPAL_COL (c) ;
1537 /* --------------------- */
1538 /* added for UMFPACK */
1539 nnewlyempty_col++ ;
1540 /* --------------------- */
1541 }
1542 else
1543 {
1544 /* set column length and set score */
1545 ASSERT (score >= 0) ;
1546 ASSERT (score <= n_col) ;
1547 Col [c].length = col_length ;
1548 Col [c].shared2.score = score ;
1549 }
1550 }
1551 DEBUG1 (("colamd: Dense, null, and newly-null columns killed: "ID"\n",
1552 n_col-n_col2)) ;
1553
1554 /* At this point, all empty rows and columns are dead. All live columns */
1555 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
1556 /* yet). Rows may contain dead columns, but all live rows contain at */
1557 /* least one live column. */
1558
1559 #ifndef NDEBUG
1560 debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
1561 #endif /* NDEBUG */
1562
1563 /* === Initialize degree lists ========================================== */
1564
1565 #ifndef NDEBUG
1566 debug_count = 0 ;
1567 #endif /* NDEBUG */
1568
1569 /* clear the hash buckets */
1570 for (c = 0 ; c <= n_col ; c++)
1571 {
1572 head [c] = EMPTY ;
1573 }
1574 min_score = n_col ;
1575 /* place in reverse order, so low column indices are at the front */
1576 /* of the lists. This is to encourage natural tie-breaking */
1577 for (c = n_col-1 ; c >= 0 ; c--)
1578 {
1579 /* only add principal columns to degree lists */
1580 if (COL_IS_ALIVE (c))
1581 {
1582 DEBUG4 (("place "ID" score "ID" minscore "ID" ncol "ID"\n",
1583 c, Col [c].shared2.score, min_score, n_col)) ;
1584
1585 /* === Add columns score to DList =============================== */
1586
1587 score = Col [c].shared2.score ;
1588
1589 ASSERT (min_score >= 0) ;
1590 ASSERT (min_score <= n_col) ;
1591 ASSERT (score >= 0) ;
1592 ASSERT (score <= n_col) ;
1593 ASSERT (head [score] >= EMPTY) ;
1594
1595 /* now add this column to dList at proper score location */
1596 next_col = head [score] ;
1597 Col [c].shared3.prev = EMPTY ;
1598 Col [c].shared4.degree_next = next_col ;
1599
1600 /* if there already was a column with the same score, set its */
1601 /* previous pointer to this new column */
1602 if (next_col != EMPTY)
1603 {
1604 Col [next_col].shared3.prev = c ;
1605 }
1606 head [score] = c ;
1607
1608 /* see if this score is less than current min */
1609 min_score = MIN (min_score, score) ;
1610
1611 #ifndef NDEBUG
1612 debug_count++ ;
1613 #endif /* NDEBUG */
1614
1615 }
1616 }
1617
1618 #ifndef NDEBUG
1619 DEBUG1 (("colamd: Live cols "ID" out of "ID", non-princ: "ID"\n",
1620 debug_count, n_col, n_col-debug_count)) ;
1621 ASSERT (debug_count == n_col2) ;
1622 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
1623 #endif /* NDEBUG */
1624
1625 /* === Return number of remaining columns, and max row degree =========== */
1626
1627 *p_n_col2 = n_col2 ;
1628 *p_n_row2 = n_row2 ;
1629 *p_max_deg = max_deg ;
1630
1631 /* --------------------- */
1632 /* added for UMFPACK */
1633 *p_ndense_row = ndense_row ;
1634 *p_nempty_row = nempty_row ; /* original empty rows */
1635 *p_nnewlyempty_row = nnewlyempty_row ;
1636 *p_ndense_col = ndense_col ;
1637 *p_nempty_col = nempty_col ; /* original empty cols */
1638 *p_nnewlyempty_col = nnewlyempty_col ;
1639 /* --------------------- */
1640 }
1641
1642
1643 /* ========================================================================== */
1644 /* === find_ordering ======================================================== */
1645 /* ========================================================================== */
1646
1647 /*
1648 Order the principal columns of the supercolumn form of the matrix
1649 (no supercolumns on input). Uses a minimum approximate column minimum
1650 degree ordering method. Not user-callable.
1651 */
1652
find_ordering(Int n_row,Int n_col,Int Alen,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int head[],Int n_col2,Int max_deg,Int pfree,Int Front_npivcol[],Int Front_nrows[],Int Front_ncols[],Int Front_parent[],Int Front_cols[],Int * p_nfr,Int aggressive,Int InFront[])1653 PRIVATE Int find_ordering /* return the number of garbage collections */
1654 (
1655 /* === Parameters ======================================================= */
1656
1657 Int n_row, /* number of rows of A */
1658 Int n_col, /* number of columns of A */
1659 Int Alen, /* size of A, 2*nnz + n_col or larger */
1660 Colamd_Row Row [], /* of size n_row+1 */
1661 Colamd_Col Col [], /* of size n_col+1 */
1662 Int A [], /* column form and row form of A */
1663 Int head [], /* of size n_col+1 */
1664 Int n_col2, /* Remaining columns to order */
1665 Int max_deg, /* Maximum row degree */
1666 Int pfree /* index of first free slot (2*nnz on entry) */
1667 /* ------------------ */
1668 /* added for UMFPACK: */
1669 , Int Front_npivcol [ ]
1670 , Int Front_nrows [ ]
1671 , Int Front_ncols [ ]
1672 , Int Front_parent [ ]
1673 , Int Front_cols [ ]
1674 , Int *p_nfr /* number of fronts */
1675 , Int aggressive
1676 , Int InFront [ ]
1677 /* ------------------ */
1678 )
1679 {
1680 /* === Local variables ================================================== */
1681
1682 Int k ; /* current pivot ordering step */
1683 Int pivot_col ; /* current pivot column */
1684 Int *cp ; /* a column pointer */
1685 Int *rp ; /* a row pointer */
1686 Int pivot_row ; /* current pivot row */
1687 Int *new_cp ; /* modified column pointer */
1688 Int *new_rp ; /* modified row pointer */
1689 Int pivot_row_start ; /* pointer to start of pivot row */
1690 Int pivot_row_degree ; /* number of columns in pivot row */
1691 Int pivot_row_length ; /* number of supercolumns in pivot row */
1692 Int pivot_col_score ; /* score of pivot column */
1693 Int needed_memory ; /* free space needed for pivot row */
1694 Int *cp_end ; /* pointer to the end of a column */
1695 Int *rp_end ; /* pointer to the end of a row */
1696 Int row ; /* a row index */
1697 Int col ; /* a column index */
1698 Int max_score ; /* maximum possible score */
1699 Int cur_score ; /* score of current column */
1700 unsigned Int hash ; /* hash value for supernode detection */
1701 Int head_column ; /* head of hash bucket */
1702 Int first_col ; /* first column in hash bucket */
1703 Int tag_mark ; /* marker value for mark array */
1704 Int row_mark ; /* Row [row].shared2.mark */
1705 Int set_difference ; /* set difference size of row with pivot row */
1706 Int min_score ; /* smallest column score */
1707 Int col_thickness ; /* "thickness" (no. of columns in a supercol) */
1708 Int max_mark ; /* maximum value of tag_mark */
1709 Int pivot_col_thickness ; /* number of columns represented by pivot col */
1710 Int prev_col ; /* Used by Dlist operations. */
1711 Int next_col ; /* Used by Dlist operations. */
1712 Int ngarbage ; /* number of garbage collections performed */
1713
1714 #ifndef NDEBUG
1715 Int debug_d ; /* debug loop counter */
1716 Int debug_step = 0 ; /* debug loop counter */
1717 #endif /* NDEBUG */
1718
1719 /* ------------------ */
1720 /* added for UMFPACK: */
1721 Int pivot_row_thickness ; /* number of rows represented by pivot row */
1722 Int nfr = 0 ; /* number of fronts */
1723 Int child ;
1724 /* ------------------ */
1725
1726 /* === Initialization and clear mark ==================================== */
1727
1728 max_mark = MAX_MARK (n_col) ; /* defined in umfpack.h */
1729 tag_mark = clear_mark (n_row, Row) ;
1730 min_score = 0 ;
1731 ngarbage = 0 ;
1732 DEBUG1 (("colamd: Ordering, n_col2="ID"\n", n_col2)) ;
1733
1734 for (row = 0 ; row < n_row ; row++)
1735 {
1736 InFront [row] = EMPTY ;
1737 }
1738
1739 /* === Order the columns ================================================ */
1740
1741 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
1742 {
1743
1744 #ifndef NDEBUG
1745 if (debug_step % 100 == 0)
1746 {
1747 DEBUG2 (("\n... Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
1748 }
1749 else
1750 {
1751 DEBUG3 (("\n-----Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
1752 }
1753 debug_step++ ;
1754 debug_deg_lists (n_row, n_col, Row, Col, head,
1755 min_score, n_col2-k, max_deg) ;
1756 debug_matrix (n_row, n_col, Row, Col, A) ;
1757 #endif /* NDEBUG */
1758
1759 /* === Select pivot column, and order it ============================ */
1760
1761 /* make sure degree list isn't empty */
1762 ASSERT (min_score >= 0) ;
1763 ASSERT (min_score <= n_col) ;
1764 ASSERT (head [min_score] >= EMPTY) ;
1765
1766 #ifndef NDEBUG
1767 for (debug_d = 0 ; debug_d < min_score ; debug_d++)
1768 {
1769 ASSERT (head [debug_d] == EMPTY) ;
1770 }
1771 #endif /* NDEBUG */
1772
1773 /* get pivot column from head of minimum degree list */
1774 while (head [min_score] == EMPTY && min_score < n_col)
1775 {
1776 min_score++ ;
1777 }
1778 pivot_col = head [min_score] ;
1779 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
1780 next_col = Col [pivot_col].shared4.degree_next ;
1781 head [min_score] = next_col ;
1782 if (next_col != EMPTY)
1783 {
1784 Col [next_col].shared3.prev = EMPTY ;
1785 }
1786
1787 ASSERT (COL_IS_ALIVE (pivot_col)) ;
1788 DEBUG3 (("Pivot col: "ID"\n", pivot_col)) ;
1789
1790 /* remember score for defrag check */
1791 pivot_col_score = Col [pivot_col].shared2.score ;
1792
1793 /* the pivot column is the kth column in the pivot order */
1794 Col [pivot_col].shared2.order = k ;
1795
1796 /* increment order count by column thickness */
1797 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
1798 /* ------------------ */
1799 /* changed for UMFPACK: */
1800 k += pivot_col_thickness ;
1801 /* ------------------ */
1802 ASSERT (pivot_col_thickness > 0) ;
1803
1804 /* === Garbage_collection, if necessary ============================= */
1805
1806 needed_memory = MIN (pivot_col_score, n_col - k) ;
1807 if (pfree + needed_memory >= Alen)
1808 {
1809 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
1810 ngarbage++ ;
1811 /* after garbage collection we will have enough */
1812 ASSERT (pfree + needed_memory < Alen) ;
1813 /* garbage collection has wiped out the Row[].shared2.mark array */
1814 tag_mark = clear_mark (n_row, Row) ;
1815
1816 #ifndef NDEBUG
1817 debug_matrix (n_row, n_col, Row, Col, A) ;
1818 #endif /* NDEBUG */
1819 }
1820
1821 /* === Compute pivot row pattern ==================================== */
1822
1823 /* get starting location for this new merged row */
1824 pivot_row_start = pfree ;
1825
1826 /* initialize new row counts to zero */
1827 pivot_row_degree = 0 ;
1828
1829 /* ------------------ */
1830 /* added for UMFPACK: */
1831 pivot_row_thickness = 0 ;
1832 /* ------------------ */
1833
1834 /* [ tag pivot column as having been visited so it isn't included */
1835 /* in merged pivot row */
1836 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
1837
1838 /* pivot row is the union of all rows in the pivot column pattern */
1839 cp = &A [Col [pivot_col].start] ;
1840 cp_end = cp + Col [pivot_col].length ;
1841 while (cp < cp_end)
1842 {
1843 /* get a row */
1844 row = *cp++ ;
1845 DEBUG4 (("Pivot col pattern %d "ID"\n", ROW_IS_ALIVE(row), row)) ;
1846 /* skip if row is dead */
1847 if (ROW_IS_DEAD (row))
1848 {
1849 continue ;
1850 }
1851
1852 /* ------------------ */
1853 /* added for UMFPACK: */
1854 /* sum the thicknesses of all the rows */
1855 /* ASSERT (Row [row].thickness > 0) ; */
1856 pivot_row_thickness += Row [row].thickness ;
1857 /* ------------------ */
1858
1859 rp = &A [Row [row].start] ;
1860 rp_end = rp + Row [row].length ;
1861 while (rp < rp_end)
1862 {
1863 /* get a column */
1864 col = *rp++ ;
1865 /* add the column, if alive and untagged */
1866 col_thickness = Col [col].shared1.thickness ;
1867 if (col_thickness > 0 && COL_IS_ALIVE (col))
1868 {
1869 /* tag column in pivot row */
1870 Col [col].shared1.thickness = -col_thickness ;
1871 ASSERT (pfree < Alen) ;
1872 /* place column in pivot row */
1873 A [pfree++] = col ;
1874 pivot_row_degree += col_thickness ;
1875 /* ------------------ */
1876 /* added for UMFPACK: */
1877 DEBUG4 (("\t\t\tNew live column in pivot row: "ID"\n",col));
1878 /* ------------------ */
1879 }
1880 /* ------------------ */
1881 /* added for UMFPACK */
1882 #ifndef NDEBUG
1883 if (col_thickness < 0 && COL_IS_ALIVE (col))
1884 {
1885 DEBUG4 (("\t\t\tOld live column in pivot row: "ID"\n",col));
1886 }
1887 #endif
1888 /* ------------------ */
1889 }
1890 }
1891
1892 /* ------------------ */
1893 /* added for UMFPACK: */
1894 /* pivot_row_thickness is the number of rows in frontal matrix */
1895 /* both pivotal rows and nonpivotal rows */
1896 /* ------------------ */
1897
1898 /* clear tag on pivot column */
1899 Col [pivot_col].shared1.thickness = pivot_col_thickness ; /* ] */
1900 max_deg = MAX (max_deg, pivot_row_degree) ;
1901
1902 #ifndef NDEBUG
1903 DEBUG3 (("check2\n")) ;
1904 debug_mark (n_row, Row, tag_mark, max_mark) ;
1905 #endif /* NDEBUG */
1906
1907 /* === Kill all rows used to construct pivot row ==================== */
1908
1909 /* also kill pivot row, temporarily */
1910 cp = &A [Col [pivot_col].start] ;
1911 cp_end = cp + Col [pivot_col].length ;
1912 while (cp < cp_end)
1913 {
1914 /* may be killing an already dead row */
1915 row = *cp++ ;
1916
1917 DEBUG2 (("Kill row in pivot col: "ID" alive? %d, front "ID"\n",
1918 row, ROW_IS_ALIVE (row), Row [row].front)) ;
1919
1920 /* added for UMFPACK: */
1921 if (ROW_IS_ALIVE (row))
1922 {
1923 if (Row [row].front != EMPTY)
1924 {
1925 /* This row represents a frontal matrix. */
1926 /* Row [row].front is a child of current front */
1927 child = Row [row].front ;
1928 Front_parent [child] = nfr ;
1929 DEBUG1 (("Front "ID" => front "ID", normal\n", child, nfr));
1930 }
1931 else
1932 {
1933 /* This is an original row. Keep track of which front
1934 * is its parent in the row-merge tree. */
1935 InFront [row] = nfr ;
1936 DEBUG1 (("Row "ID" => front "ID", normal\n", row, nfr)) ;
1937 }
1938 }
1939
1940 KILL_ROW (row) ;
1941
1942 /* ------------------ */
1943 /* added for UMFPACK: */
1944 Row [row].thickness = 0 ;
1945 /* ------------------ */
1946 }
1947
1948 /* === Select a row index to use as the new pivot row =============== */
1949
1950 pivot_row_length = pfree - pivot_row_start ;
1951 if (pivot_row_length > 0)
1952 {
1953 /* pick the "pivot" row arbitrarily (first row in col) */
1954 pivot_row = A [Col [pivot_col].start] ;
1955 DEBUG3 (("Pivotal row is "ID"\n", pivot_row)) ;
1956 }
1957 else
1958 {
1959 /* there is no pivot row, since it is of zero length */
1960 pivot_row = EMPTY ;
1961 ASSERT (pivot_row_length == 0) ;
1962 }
1963 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
1964
1965 /* === Approximate degree computation =============================== */
1966
1967 /* Here begins the computation of the approximate degree. The column */
1968 /* score is the sum of the pivot row "length", plus the size of the */
1969 /* set differences of each row in the column minus the pattern of the */
1970 /* pivot row itself. The column ("thickness") itself is also */
1971 /* excluded from the column score (we thus use an approximate */
1972 /* external degree). */
1973
1974 /* The time taken by the following code (compute set differences, and */
1975 /* add them up) is proportional to the size of the data structure */
1976 /* being scanned - that is, the sum of the sizes of each column in */
1977 /* the pivot row. Thus, the amortized time to compute a column score */
1978 /* is proportional to the size of that column (where size, in this */
1979 /* context, is the column "length", or the number of row indices */
1980 /* in that column). The number of row indices in a column is */
1981 /* monotonically non-decreasing, from the length of the original */
1982 /* column on input to colamd. */
1983
1984 /* === Compute set differences ====================================== */
1985
1986 DEBUG3 (("** Computing set differences phase. **\n")) ;
1987
1988 /* pivot row is currently dead - it will be revived later. */
1989
1990 DEBUG3 (("Pivot row: \n")) ;
1991 /* for each column in pivot row */
1992 rp = &A [pivot_row_start] ;
1993 rp_end = rp + pivot_row_length ;
1994 while (rp < rp_end)
1995 {
1996 col = *rp++ ;
1997 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1998 DEBUG3 ((" Col: "ID"\n", col)) ;
1999
2000 /* clear tags used to construct pivot row pattern */
2001 col_thickness = -Col [col].shared1.thickness ;
2002 ASSERT (col_thickness > 0) ;
2003 Col [col].shared1.thickness = col_thickness ;
2004
2005 /* === Remove column from degree list =========================== */
2006
2007 cur_score = Col [col].shared2.score ;
2008 prev_col = Col [col].shared3.prev ;
2009 next_col = Col [col].shared4.degree_next ;
2010 ASSERT (cur_score >= 0) ;
2011 ASSERT (cur_score <= n_col) ;
2012 ASSERT (cur_score >= EMPTY) ;
2013 if (prev_col == EMPTY)
2014 {
2015 head [cur_score] = next_col ;
2016 }
2017 else
2018 {
2019 Col [prev_col].shared4.degree_next = next_col ;
2020 }
2021 if (next_col != EMPTY)
2022 {
2023 Col [next_col].shared3.prev = prev_col ;
2024 }
2025
2026 /* === Scan the column ========================================== */
2027
2028 cp = &A [Col [col].start] ;
2029 cp_end = cp + Col [col].length ;
2030 while (cp < cp_end)
2031 {
2032 /* get a row */
2033 row = *cp++ ;
2034 row_mark = Row [row].shared2.mark ;
2035 /* skip if dead */
2036 if (ROW_IS_MARKED_DEAD (row_mark))
2037 {
2038 continue ;
2039 }
2040 ASSERT (row != pivot_row) ;
2041 set_difference = row_mark - tag_mark ;
2042 /* check if the row has been seen yet */
2043 if (set_difference < 0)
2044 {
2045 ASSERT (Row [row].shared1.degree <= max_deg) ;
2046 set_difference = Row [row].shared1.degree ;
2047 }
2048 /* subtract column thickness from this row's set difference */
2049 set_difference -= col_thickness ;
2050 ASSERT (set_difference >= 0) ;
2051 ASSERT (ROW_IS_ALIVE (row)) ;
2052
2053 /* absorb this row if the set difference becomes zero */
2054 if (set_difference == 0 && aggressive)
2055 {
2056 /* v4.1: do aggressive absorption */
2057 DEBUG3 (("aggressive absorption. Row: "ID"\n", row)) ;
2058
2059 if (Row [row].front != EMPTY)
2060 {
2061 /* Row [row].front is a child of current front. */
2062 child = Row [row].front ;
2063 Front_parent [child] = nfr ;
2064 DEBUG1 (("Front "ID" => front "ID", aggressive\n",
2065 child, nfr)) ;
2066 }
2067 else
2068 {
2069 /* this is an original row. Keep track of which front
2070 * assembles it, for the row-merge tree */
2071 InFront [row] = nfr ;
2072 DEBUG1 (("Row "ID" => front "ID", aggressive\n",
2073 row, nfr)) ;
2074 }
2075
2076 KILL_ROW (row) ;
2077
2078 /* sum the thicknesses of all the rows */
2079 /* ASSERT (Row [row].thickness > 0) ; */
2080 pivot_row_thickness += Row [row].thickness ;
2081 Row [row].thickness = 0 ;
2082
2083 }
2084 else
2085 {
2086 /* save the new mark */
2087 Row [row].shared2.mark = set_difference + tag_mark ;
2088 }
2089 }
2090 }
2091
2092 #ifndef NDEBUG
2093 debug_deg_lists (n_row, n_col, Row, Col, head,
2094 min_score, n_col2-k-pivot_row_degree, max_deg) ;
2095 #endif /* NDEBUG */
2096
2097 /* === Add up set differences for each column ======================= */
2098
2099 DEBUG3 (("** Adding set differences phase. **\n")) ;
2100
2101 /* for each column in pivot row */
2102 rp = &A [pivot_row_start] ;
2103 rp_end = rp + pivot_row_length ;
2104 while (rp < rp_end)
2105 {
2106 /* get a column */
2107 col = *rp++ ;
2108 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2109 hash = 0 ;
2110 cur_score = 0 ;
2111 cp = &A [Col [col].start] ;
2112 /* compact the column */
2113 new_cp = cp ;
2114 cp_end = cp + Col [col].length ;
2115
2116 DEBUG4 (("Adding set diffs for Col: "ID".\n", col)) ;
2117
2118 while (cp < cp_end)
2119 {
2120 /* get a row */
2121 row = *cp++ ;
2122 ASSERT(row >= 0 && row < n_row) ;
2123 row_mark = Row [row].shared2.mark ;
2124 /* skip if dead */
2125 if (ROW_IS_MARKED_DEAD (row_mark))
2126 {
2127 /* ------------------ */
2128 /* changed for UMFPACK: */
2129 DEBUG4 ((" Row "ID", dead\n", row)) ;
2130 /* ------------------ */
2131 continue ;
2132 }
2133 /* ------------------ */
2134 /* changed for UMFPACK: */
2135 /* ASSERT (row_mark > tag_mark) ; */
2136 DEBUG4 ((" Row "ID", set diff "ID"\n", row, row_mark-tag_mark));
2137 ASSERT (row_mark >= tag_mark) ;
2138 /* ------------------ */
2139 /* compact the column */
2140 *new_cp++ = row ;
2141 /* compute hash function */
2142 hash += row ;
2143 /* add set difference */
2144 cur_score += row_mark - tag_mark ;
2145 /* integer overflow... */
2146 cur_score = MIN (cur_score, n_col) ;
2147 }
2148
2149 /* recompute the column's length */
2150 Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
2151
2152 /* === Further mass elimination ================================= */
2153
2154 if (Col [col].length == 0)
2155 {
2156 DEBUG4 (("further mass elimination. Col: "ID"\n", col)) ;
2157 /* nothing left but the pivot row in this column */
2158 KILL_PRINCIPAL_COL (col) ;
2159 pivot_row_degree -= Col [col].shared1.thickness ;
2160 ASSERT (pivot_row_degree >= 0) ;
2161 /* order it */
2162 Col [col].shared2.order = k ;
2163 /* increment order count by column thickness */
2164 k += Col [col].shared1.thickness ;
2165
2166 /* ------------------ */
2167 /* added for UMFPACK: */
2168 pivot_col_thickness += Col [col].shared1.thickness ;
2169
2170 /* add to column list of front ... */
2171 #ifndef NDEBUG
2172 DEBUG1 (("Mass")) ;
2173 dump_super (col, Col, n_col) ;
2174 #endif
2175 Col [Col [col].lastcol].nextcol = Front_cols [nfr] ;
2176 Front_cols [nfr] = col ;
2177 /* ------------------ */
2178
2179 }
2180 else
2181 {
2182 /* === Prepare for supercolumn detection ==================== */
2183
2184 DEBUG4 (("Preparing supercol detection for Col: "ID".\n", col));
2185
2186 /* save score so far */
2187 Col [col].shared2.score = cur_score ;
2188
2189 /* add column to hash table, for supercolumn detection */
2190 /* NOTE: hash is an unsigned Int to avoid a problem in ANSI C.
2191 * The sign of the expression a % b is not defined when a and/or
2192 * b are negative. Since hash is unsigned and n_col >= 0,
2193 * this problem is avoided. */
2194 hash %= n_col + 1 ;
2195
2196 DEBUG4 ((" Hash = "ID", n_col = "ID".\n", (Int) hash, n_col)) ;
2197 ASSERT (((Int) hash) <= n_col) ;
2198
2199 head_column = head [hash] ;
2200 if (head_column > EMPTY)
2201 {
2202 /* degree list "hash" is non-empty, use prev (shared3) of */
2203 /* first column in degree list as head of hash bucket */
2204 first_col = Col [head_column].shared3.headhash ;
2205 Col [head_column].shared3.headhash = col ;
2206 }
2207 else
2208 {
2209 /* degree list "hash" is empty, use head as hash bucket */
2210 first_col = - (head_column + 2) ;
2211 head [hash] = - (col + 2) ;
2212 }
2213 Col [col].shared4.hash_next = first_col ;
2214
2215 /* save hash function in Col [col].shared3.hash */
2216 Col [col].shared3.hash = (Int) hash ;
2217 ASSERT (COL_IS_ALIVE (col)) ;
2218 }
2219 }
2220
2221 /* The approximate external column degree is now computed. */
2222
2223 /* === Supercolumn detection ======================================== */
2224
2225 DEBUG3 (("** Supercolumn detection phase. **\n")) ;
2226
2227 detect_super_cols (
2228
2229 #ifndef NDEBUG
2230 n_col, Row,
2231 #endif /* NDEBUG */
2232
2233 Col, A, head, pivot_row_start, pivot_row_length) ;
2234
2235 /* === Kill the pivotal column ====================================== */
2236
2237 KILL_PRINCIPAL_COL (pivot_col) ;
2238
2239 /* ------------------ */
2240 /* added for UMFPACK: */
2241 /* add columns to column list of front */
2242 #ifndef NDEBUG
2243 DEBUG1 (("Pivot")) ;
2244 dump_super (pivot_col, Col, n_col) ;
2245 #endif
2246 Col [Col [pivot_col].lastcol].nextcol = Front_cols [nfr] ;
2247 Front_cols [nfr] = pivot_col ;
2248 /* ------------------ */
2249
2250 /* === Clear mark =================================================== */
2251
2252 tag_mark += (max_deg + 1) ;
2253 if (tag_mark >= max_mark)
2254 {
2255 DEBUG2 (("clearing tag_mark\n")) ;
2256 tag_mark = clear_mark (n_row, Row) ;
2257 }
2258
2259 #ifndef NDEBUG
2260 DEBUG3 (("check3\n")) ;
2261 debug_mark (n_row, Row, tag_mark, max_mark) ;
2262 #endif /* NDEBUG */
2263
2264 /* === Finalize the new pivot row, and column scores ================ */
2265
2266 DEBUG3 (("** Finalize scores phase. **\n")) ;
2267 DEBUG3 (("pivot_row_degree "ID"\n", pivot_row_degree)) ;
2268
2269 /* for each column in pivot row */
2270 rp = &A [pivot_row_start] ;
2271 /* compact the pivot row */
2272 new_rp = rp ;
2273 rp_end = rp + pivot_row_length ;
2274 while (rp < rp_end)
2275 {
2276 col = *rp++ ;
2277 DEBUG3 (("Col "ID" \n", col)) ;
2278 /* skip dead columns */
2279 if (COL_IS_DEAD (col))
2280 {
2281 DEBUG3 (("dead\n")) ;
2282 continue ;
2283 }
2284 *new_rp++ = col ;
2285 /* add new pivot row to column */
2286 A [Col [col].start + (Col [col].length++)] = pivot_row ;
2287
2288 /* retrieve score so far and add on pivot row's degree. */
2289 /* (we wait until here for this in case the pivot */
2290 /* row's degree was reduced due to mass elimination). */
2291 cur_score = Col [col].shared2.score + pivot_row_degree ;
2292 DEBUG3 ((" cur_score "ID" ", cur_score)) ;
2293
2294 /* calculate the max possible score as the number of */
2295 /* external columns minus the 'k' value minus the */
2296 /* columns thickness */
2297 max_score = n_col - k - Col [col].shared1.thickness ;
2298 DEBUG3 ((" max_score "ID" ", max_score)) ;
2299
2300 /* make the score the external degree of the union-of-rows */
2301 cur_score -= Col [col].shared1.thickness ;
2302 DEBUG3 ((" cur_score "ID" ", cur_score)) ;
2303
2304 /* make sure score is less or equal than the max score */
2305 cur_score = MIN (cur_score, max_score) ;
2306 ASSERT (cur_score >= 0) ;
2307
2308 /* store updated score */
2309 Col [col].shared2.score = cur_score ;
2310 DEBUG3 ((" "ID"\n", cur_score)) ;
2311
2312 /* === Place column back in degree list ========================= */
2313
2314 ASSERT (min_score >= 0) ;
2315 ASSERT (min_score <= n_col) ;
2316 ASSERT (cur_score >= 0) ;
2317 ASSERT (cur_score <= n_col) ;
2318 ASSERT (head [cur_score] >= EMPTY) ;
2319 next_col = head [cur_score] ;
2320 Col [col].shared4.degree_next = next_col ;
2321 Col [col].shared3.prev = EMPTY ;
2322 if (next_col != EMPTY)
2323 {
2324 Col [next_col].shared3.prev = col ;
2325 }
2326 head [cur_score] = col ;
2327
2328 /* see if this score is less than current min */
2329 min_score = MIN (min_score, cur_score) ;
2330
2331 }
2332
2333 #ifndef NDEBUG
2334 debug_deg_lists (n_row, n_col, Row, Col, head,
2335 min_score, n_col2-k, max_deg) ;
2336 #endif /* NDEBUG */
2337
2338 /* ------------------ */
2339 /* added for UMFPACK: */
2340 /* frontal matrix can have more pivot cols than pivot rows for */
2341 /* singular matrices. */
2342
2343 /* number of candidate pivot columns */
2344 Front_npivcol [nfr] = pivot_col_thickness ;
2345
2346 /* all rows (not just size of contrib. block) */
2347 Front_nrows [nfr] = pivot_row_thickness ;
2348
2349 /* all cols */
2350 Front_ncols [nfr] = pivot_col_thickness + pivot_row_degree ;
2351
2352 Front_parent [nfr] = EMPTY ;
2353
2354 pivot_row_thickness -= pivot_col_thickness ;
2355 DEBUG1 (("Front "ID" Pivot_row_thickness after pivot cols elim: "ID"\n",
2356 nfr, pivot_row_thickness)) ;
2357 pivot_row_thickness = MAX (0, pivot_row_thickness) ;
2358 /* ------------------ */
2359
2360 /* === Resurrect the new pivot row ================================== */
2361
2362 if (pivot_row_degree > 0
2363 /* ------------------ */
2364 /* added for UMFPACK. Note that this part of the expression should be
2365 * removed if this routine is used outside of UMFPACK, for a Cholesky
2366 * factorization of (AQ)'(AQ) */
2367 && pivot_row_thickness > 0
2368 /* ------------------ */
2369 )
2370 {
2371 /* update pivot row length to reflect any cols that were killed */
2372 /* during super-col detection and mass elimination */
2373 Row [pivot_row].start = pivot_row_start ;
2374 Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
2375 ASSERT (Row [pivot_row].length > 0) ;
2376 Row [pivot_row].shared1.degree = pivot_row_degree ;
2377 Row [pivot_row].shared2.mark = 0 ;
2378 /* ------------------ */
2379 /* added for UMFPACK: */
2380 Row [pivot_row].thickness = pivot_row_thickness ;
2381 Row [pivot_row].front = nfr ;
2382 /* ------------------ */
2383 /* pivot row is no longer dead */
2384 }
2385
2386 /* ------------------ */
2387 /* added for UMFPACK: */
2388
2389 #ifndef NDEBUG
2390 DEBUG1 (("Front "ID" : "ID" "ID" "ID" ", nfr,
2391 Front_npivcol [nfr], Front_nrows [nfr], Front_ncols [nfr])) ;
2392 DEBUG1 ((" cols:[ ")) ;
2393 debug_d = 0 ;
2394 for (col = Front_cols [nfr] ; col != EMPTY ; col = Col [col].nextcol)
2395 {
2396 DEBUG1 ((" "ID, col)) ;
2397 ASSERT (col >= 0 && col < n_col) ;
2398 ASSERT (COL_IS_DEAD (col)) ;
2399 debug_d++ ;
2400 ASSERT (debug_d <= pivot_col_thickness) ;
2401 }
2402 ASSERT (debug_d == pivot_col_thickness) ;
2403 DEBUG1 ((" ]\n ")) ;
2404 #endif
2405 nfr++ ; /* one more front */
2406 /* ------------------ */
2407
2408 }
2409
2410 /* === All principal columns have now been ordered ====================== */
2411
2412 /* ------------------ */
2413 /* added for UMFPACK: */
2414 *p_nfr = nfr ;
2415 /* ------------------ */
2416
2417 return (ngarbage) ;
2418 }
2419
2420
2421 /* ========================================================================== */
2422 /* === order_children deleted for UMFPACK =================================== */
2423 /* ========================================================================== */
2424
2425 /* ========================================================================== */
2426 /* === detect_super_cols ==================================================== */
2427 /* ========================================================================== */
2428
2429 /*
2430 Detects supercolumns by finding matches between columns in the hash buckets.
2431 Check amongst columns in the set A [row_start ... row_start + row_length-1].
2432 The columns under consideration are currently *not* in the degree lists,
2433 and have already been placed in the hash buckets.
2434
2435 The hash bucket for columns whose hash function is equal to h is stored
2436 as follows:
2437
2438 if head [h] is >= 0, then head [h] contains a degree list, so:
2439
2440 head [h] is the first column in degree bucket h.
2441 Col [head [h]].headhash gives the first column in hash bucket h.
2442
2443 otherwise, the degree list is empty, and:
2444
2445 -(head [h] + 2) is the first column in hash bucket h.
2446
2447 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
2448 column" pointer. Col [c].shared3.hash is used instead as the hash number
2449 for that column. The value of Col [c].shared4.hash_next is the next column
2450 in the same hash bucket.
2451
2452 Assuming no, or "few" hash collisions, the time taken by this routine is
2453 linear in the sum of the sizes (lengths) of each column whose score has
2454 just been computed in the approximate degree computation.
2455 Not user-callable.
2456 */
2457
detect_super_cols(Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int head[],Int row_start,Int row_length)2458 PRIVATE void detect_super_cols
2459 (
2460 /* === Parameters ======================================================= */
2461
2462 #ifndef NDEBUG
2463 /* these two parameters are only needed when debugging is enabled: */
2464 Int n_col, /* number of columns of A */
2465 Colamd_Row Row [], /* of size n_row+1 */
2466 #endif /* NDEBUG */
2467
2468 Colamd_Col Col [], /* of size n_col+1 */
2469 Int A [], /* row indices of A */
2470 Int head [], /* head of degree lists and hash buckets */
2471 Int row_start, /* pointer to set of columns to check */
2472 Int row_length /* number of columns to check */
2473 )
2474 {
2475 /* === Local variables ================================================== */
2476
2477 Int hash ; /* hash value for a column */
2478 Int *rp ; /* pointer to a row */
2479 Int c ; /* a column index */
2480 Int super_c ; /* column index of the column to absorb into */
2481 Int *cp1 ; /* column pointer for column super_c */
2482 Int *cp2 ; /* column pointer for column c */
2483 Int length ; /* length of column super_c */
2484 Int prev_c ; /* column preceding c in hash bucket */
2485 Int i ; /* loop counter */
2486 Int *rp_end ; /* pointer to the end of the row */
2487 Int col ; /* a column index in the row to check */
2488 Int head_column ; /* first column in hash bucket or degree list */
2489 Int first_col ; /* first column in hash bucket */
2490
2491 /* === Consider each column in the row ================================== */
2492
2493 rp = &A [row_start] ;
2494 rp_end = rp + row_length ;
2495 while (rp < rp_end)
2496 {
2497 col = *rp++ ;
2498 if (COL_IS_DEAD (col))
2499 {
2500 continue ;
2501 }
2502
2503 /* get hash number for this column */
2504 hash = Col [col].shared3.hash ;
2505 ASSERT (hash <= n_col) ;
2506
2507 /* === Get the first column in this hash bucket ===================== */
2508
2509 head_column = head [hash] ;
2510 if (head_column > EMPTY)
2511 {
2512 first_col = Col [head_column].shared3.headhash ;
2513 }
2514 else
2515 {
2516 first_col = - (head_column + 2) ;
2517 }
2518
2519 /* === Consider each column in the hash bucket ====================== */
2520
2521 for (super_c = first_col ; super_c != EMPTY ;
2522 super_c = Col [super_c].shared4.hash_next)
2523 {
2524 ASSERT (COL_IS_ALIVE (super_c)) ;
2525 ASSERT (Col [super_c].shared3.hash == hash) ;
2526 length = Col [super_c].length ;
2527
2528 /* prev_c is the column preceding column c in the hash bucket */
2529 prev_c = super_c ;
2530
2531 /* === Compare super_c with all columns after it ================ */
2532
2533 for (c = Col [super_c].shared4.hash_next ;
2534 c != EMPTY ; c = Col [c].shared4.hash_next)
2535 {
2536 ASSERT (c != super_c) ;
2537 ASSERT (COL_IS_ALIVE (c)) ;
2538 ASSERT (Col [c].shared3.hash == hash) ;
2539
2540 /* not identical if lengths or scores are different */
2541 if (Col [c].length != length ||
2542 Col [c].shared2.score != Col [super_c].shared2.score)
2543 {
2544 prev_c = c ;
2545 continue ;
2546 }
2547
2548 /* compare the two columns */
2549 cp1 = &A [Col [super_c].start] ;
2550 cp2 = &A [Col [c].start] ;
2551
2552 for (i = 0 ; i < length ; i++)
2553 {
2554 /* the columns are "clean" (no dead rows) */
2555 ASSERT (ROW_IS_ALIVE (*cp1)) ;
2556 ASSERT (ROW_IS_ALIVE (*cp2)) ;
2557 /* row indices will same order for both supercols, */
2558 /* no gather scatter nessasary */
2559 if (*cp1++ != *cp2++)
2560 {
2561 break ;
2562 }
2563 }
2564
2565 /* the two columns are different if the for-loop "broke" */
2566 if (i != length)
2567 {
2568 prev_c = c ;
2569 continue ;
2570 }
2571
2572 /* === Got it! two columns are identical =================== */
2573
2574 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
2575
2576 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
2577 Col [c].shared1.parent = super_c ;
2578 KILL_NON_PRINCIPAL_COL (c) ;
2579
2580 Col [c].shared2.order = EMPTY ;
2581 /* remove c from hash bucket */
2582 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
2583
2584 /* ------------------ */
2585 /* added for UMFPACK: */
2586 /* add c to end of list of super_c */
2587 ASSERT (Col [super_c].lastcol >= 0) ;
2588 ASSERT (Col [super_c].lastcol < n_col) ;
2589 Col [Col [super_c].lastcol].nextcol = c ;
2590 Col [super_c].lastcol = Col [c].lastcol ;
2591 #ifndef NDEBUG
2592 /* dump the supercolumn */
2593 DEBUG1 (("Super")) ;
2594 dump_super (super_c, Col, n_col) ;
2595 #endif
2596 /* ------------------ */
2597
2598 }
2599 }
2600
2601 /* === Empty this hash bucket ======================================= */
2602
2603 if (head_column > EMPTY)
2604 {
2605 /* corresponding degree list "hash" is not empty */
2606 Col [head_column].shared3.headhash = EMPTY ;
2607 }
2608 else
2609 {
2610 /* corresponding degree list "hash" is empty */
2611 head [hash] = EMPTY ;
2612 }
2613 }
2614 }
2615
2616
2617 /* ========================================================================== */
2618 /* === garbage_collection =================================================== */
2619 /* ========================================================================== */
2620
2621 /*
2622 Defragments and compacts columns and rows in the workspace A. Used when
2623 all avaliable memory has been used while performing row merging. Returns
2624 the index of the first free position in A, after garbage collection. The
2625 time taken by this routine is linear is the size of the array A, which is
2626 itself linear in the number of nonzeros in the input matrix.
2627 Not user-callable.
2628 */
2629
garbage_collection(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int * pfree)2630 PRIVATE Int garbage_collection /* returns the new value of pfree */
2631 (
2632 /* === Parameters ======================================================= */
2633
2634 Int n_row, /* number of rows */
2635 Int n_col, /* number of columns */
2636 Colamd_Row Row [], /* row info */
2637 Colamd_Col Col [], /* column info */
2638 Int A [], /* A [0 ... Alen-1] holds the matrix */
2639 Int *pfree /* &A [0] ... pfree is in use */
2640 )
2641 {
2642 /* === Local variables ================================================== */
2643
2644 Int *psrc ; /* source pointer */
2645 Int *pdest ; /* destination pointer */
2646 Int j ; /* counter */
2647 Int r ; /* a row index */
2648 Int c ; /* a column index */
2649 Int length ; /* length of a row or column */
2650
2651 #ifndef NDEBUG
2652 Int debug_rows ;
2653 DEBUG2 (("Defrag..\n")) ;
2654 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
2655 debug_rows = 0 ;
2656 #endif /* NDEBUG */
2657
2658 /* === Defragment the columns =========================================== */
2659
2660 pdest = &A[0] ;
2661 for (c = 0 ; c < n_col ; c++)
2662 {
2663 if (COL_IS_ALIVE (c))
2664 {
2665 psrc = &A [Col [c].start] ;
2666
2667 /* move and compact the column */
2668 ASSERT (pdest <= psrc) ;
2669 Col [c].start = (Int) (pdest - &A [0]) ;
2670 length = Col [c].length ;
2671 for (j = 0 ; j < length ; j++)
2672 {
2673 r = *psrc++ ;
2674 if (ROW_IS_ALIVE (r))
2675 {
2676 *pdest++ = r ;
2677 }
2678 }
2679 Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
2680 }
2681 }
2682
2683 /* === Prepare to defragment the rows =================================== */
2684
2685 for (r = 0 ; r < n_row ; r++)
2686 {
2687 if (ROW_IS_ALIVE (r))
2688 {
2689 if (Row [r].length == 0)
2690 {
2691 /* :: defrag row kill :: */
2692 /* This row is of zero length. cannot compact it, so kill it.
2693 * NOTE: in the current version, there are no zero-length live
2694 * rows when garbage_collection is called. So this code will
2695 * never trigger. However, if the code is modified, or if
2696 * garbage_collection is called at a different place, then rows
2697 * can be of zero length. So this test is kept, just in case.
2698 */
2699 DEBUGm4 (("Defrag row kill\n")) ;
2700 KILL_ROW (r) ;
2701 }
2702 else
2703 {
2704 /* save first column index in Row [r].shared2.first_column */
2705 psrc = &A [Row [r].start] ;
2706 Row [r].shared2.first_column = *psrc ;
2707 ASSERT (ROW_IS_ALIVE (r)) ;
2708 /* flag the start of the row with the one's complement of row */
2709 *psrc = ONES_COMPLEMENT (r) ;
2710 #ifndef NDEBUG
2711 debug_rows++ ;
2712 #endif /* NDEBUG */
2713 }
2714 }
2715 }
2716
2717 /* === Defragment the rows ============================================== */
2718
2719 psrc = pdest ;
2720 while (psrc < pfree)
2721 {
2722 /* find a negative number ... the start of a row */
2723 if (*psrc++ < 0)
2724 {
2725 psrc-- ;
2726 /* get the row index */
2727 r = ONES_COMPLEMENT (*psrc) ;
2728 ASSERT (r >= 0 && r < n_row) ;
2729 /* restore first column index */
2730 *psrc = Row [r].shared2.first_column ;
2731 ASSERT (ROW_IS_ALIVE (r)) ;
2732
2733 /* move and compact the row */
2734 ASSERT (pdest <= psrc) ;
2735 Row [r].start = (Int) (pdest - &A [0]) ;
2736 length = Row [r].length ;
2737 for (j = 0 ; j < length ; j++)
2738 {
2739 c = *psrc++ ;
2740 if (COL_IS_ALIVE (c))
2741 {
2742 *pdest++ = c ;
2743 }
2744 }
2745 Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
2746
2747 #ifndef NDEBUG
2748 debug_rows-- ;
2749 #endif /* NDEBUG */
2750
2751 }
2752 }
2753 /* ensure we found all the rows */
2754 ASSERT (debug_rows == 0) ;
2755
2756 /* === Return the new value of pfree ==================================== */
2757
2758 return ((Int) (pdest - &A [0])) ;
2759 }
2760
2761
2762 /* ========================================================================== */
2763 /* === clear_mark =========================================================== */
2764 /* ========================================================================== */
2765
2766 /*
2767 Clears the Row [].shared2.mark array, and returns the new tag_mark.
2768 Return value is the new tag_mark. Not user-callable.
2769 */
2770
clear_mark(Int n_row,Colamd_Row Row[])2771 PRIVATE Int clear_mark /* return the new value for tag_mark */
2772 (
2773 /* === Parameters ======================================================= */
2774
2775 Int n_row, /* number of rows in A */
2776 Colamd_Row Row [] /* Row [0 ... n-1].shared2.mark is set to zero */
2777 )
2778 {
2779 /* === Local variables ================================================== */
2780
2781 Int r ;
2782
2783 for (r = 0 ; r < n_row ; r++)
2784 {
2785 if (ROW_IS_ALIVE (r))
2786 {
2787 Row [r].shared2.mark = 0 ;
2788 }
2789 }
2790
2791 /* ------------------ */
2792 return (1) ;
2793 /* ------------------ */
2794
2795 }
2796
2797
2798 /* ========================================================================== */
2799 /* === print_report removed for UMFPACK ===================================== */
2800 /* ========================================================================== */
2801
2802
2803
2804 /* ========================================================================== */
2805 /* === colamd debugging routines ============================================ */
2806 /* ========================================================================== */
2807
2808 /* When debugging is disabled, the remainder of this file is ignored. */
2809
2810 #ifndef NDEBUG
2811
2812
2813 /* ========================================================================== */
2814 /* === debug_structures ===================================================== */
2815 /* ========================================================================== */
2816
2817 /*
2818 At this point, all empty rows and columns are dead. All live columns
2819 are "clean" (containing no dead rows) and simplicial (no supercolumns
2820 yet). Rows may contain dead columns, but all live rows contain at
2821 least one live column.
2822 */
2823
debug_structures(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[],Int n_col2)2824 PRIVATE void debug_structures
2825 (
2826 /* === Parameters ======================================================= */
2827
2828 Int n_row,
2829 Int n_col,
2830 Colamd_Row Row [],
2831 Colamd_Col Col [],
2832 Int A [],
2833 Int n_col2
2834 )
2835 {
2836 /* === Local variables ================================================== */
2837
2838 Int i ;
2839 Int c ;
2840 Int *cp ;
2841 Int *cp_end ;
2842 Int len ;
2843 Int score ;
2844 Int r ;
2845 Int *rp ;
2846 Int *rp_end ;
2847 Int deg ;
2848
2849 /* === Check A, Row, and Col ============================================ */
2850
2851 for (c = 0 ; c < n_col ; c++)
2852 {
2853 if (COL_IS_ALIVE (c))
2854 {
2855 len = Col [c].length ;
2856 score = Col [c].shared2.score ;
2857 DEBUG4 (("initial live col "ID" "ID" "ID"\n", c, len, score)) ;
2858 ASSERT (len > 0) ;
2859 ASSERT (score >= 0) ;
2860 ASSERT (Col [c].shared1.thickness == 1) ;
2861 cp = &A [Col [c].start] ;
2862 cp_end = cp + len ;
2863 while (cp < cp_end)
2864 {
2865 r = *cp++ ;
2866 ASSERT (ROW_IS_ALIVE (r)) ;
2867 }
2868 }
2869 else
2870 {
2871 i = Col [c].shared2.order ;
2872 ASSERT (i >= n_col2 && i < n_col) ;
2873 }
2874 }
2875
2876 for (r = 0 ; r < n_row ; r++)
2877 {
2878 if (ROW_IS_ALIVE (r))
2879 {
2880 i = 0 ;
2881 len = Row [r].length ;
2882 deg = Row [r].shared1.degree ;
2883 ASSERT (len > 0) ;
2884 ASSERT (deg > 0) ;
2885 rp = &A [Row [r].start] ;
2886 rp_end = rp + len ;
2887 while (rp < rp_end)
2888 {
2889 c = *rp++ ;
2890 if (COL_IS_ALIVE (c))
2891 {
2892 i++ ;
2893 }
2894 }
2895 ASSERT (i > 0) ;
2896 }
2897 }
2898 }
2899
2900
2901 /* ========================================================================== */
2902 /* === debug_deg_lists ====================================================== */
2903 /* ========================================================================== */
2904
2905 /*
2906 Prints the contents of the degree lists. Counts the number of columns
2907 in the degree list and compares it to the total it should have. Also
2908 checks the row degrees.
2909 */
2910
debug_deg_lists(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int head[],Int min_score,Int should,Int max_deg)2911 PRIVATE void debug_deg_lists
2912 (
2913 /* === Parameters ======================================================= */
2914
2915 Int n_row,
2916 Int n_col,
2917 Colamd_Row Row [],
2918 Colamd_Col Col [],
2919 Int head [],
2920 Int min_score,
2921 Int should,
2922 Int max_deg
2923 )
2924 {
2925 /* === Local variables ================================================== */
2926
2927 Int deg ;
2928 Int col ;
2929 Int have ;
2930 Int row ;
2931
2932 /* === Check the degree lists =========================================== */
2933
2934 if (n_col > 10000 && UMF_debug <= 0)
2935 {
2936 return ;
2937 }
2938 have = 0 ;
2939 DEBUG4 (("Degree lists: "ID"\n", min_score)) ;
2940 for (deg = 0 ; deg <= n_col ; deg++)
2941 {
2942 col = head [deg] ;
2943 if (col == EMPTY)
2944 {
2945 continue ;
2946 }
2947 DEBUG4 ((ID":", deg)) ;
2948 while (col != EMPTY)
2949 {
2950 DEBUG4 ((" "ID, col)) ;
2951 have += Col [col].shared1.thickness ;
2952 ASSERT (COL_IS_ALIVE (col)) ;
2953 col = Col [col].shared4.degree_next ;
2954 }
2955 DEBUG4 (("\n")) ;
2956 }
2957 DEBUG4 (("should "ID" have "ID"\n", should, have)) ;
2958 ASSERT (should == have) ;
2959
2960 /* === Check the row degrees ============================================ */
2961
2962 if (n_row > 10000 && UMF_debug <= 0)
2963 {
2964 return ;
2965 }
2966 for (row = 0 ; row < n_row ; row++)
2967 {
2968 if (ROW_IS_ALIVE (row))
2969 {
2970 ASSERT (Row [row].shared1.degree <= max_deg) ;
2971 }
2972 }
2973 }
2974
2975
2976 /* ========================================================================== */
2977 /* === debug_mark =========================================================== */
2978 /* ========================================================================== */
2979
2980 /*
2981 Ensures that the tag_mark is less that the maximum and also ensures that
2982 each entry in the mark array is less than the tag mark.
2983 */
2984
debug_mark(Int n_row,Colamd_Row Row[],Int tag_mark,Int max_mark)2985 PRIVATE void debug_mark
2986 (
2987 /* === Parameters ======================================================= */
2988
2989 Int n_row,
2990 Colamd_Row Row [],
2991 Int tag_mark,
2992 Int max_mark
2993 )
2994 {
2995 /* === Local variables ================================================== */
2996
2997 Int r ;
2998
2999 /* === Check the Row marks ============================================== */
3000
3001 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
3002 if (n_row > 10000 && UMF_debug <= 0)
3003 {
3004 return ;
3005 }
3006 for (r = 0 ; r < n_row ; r++)
3007 {
3008 ASSERT (Row [r].shared2.mark < tag_mark) ;
3009 }
3010 }
3011
3012
3013 /* ========================================================================== */
3014 /* === debug_matrix ========================================================= */
3015 /* ========================================================================== */
3016
3017 /*
3018 Prints out the contents of the columns and the rows.
3019 */
3020
debug_matrix(Int n_row,Int n_col,Colamd_Row Row[],Colamd_Col Col[],Int A[])3021 PRIVATE void debug_matrix
3022 (
3023 /* === Parameters ======================================================= */
3024
3025 Int n_row,
3026 Int n_col,
3027 Colamd_Row Row [],
3028 Colamd_Col Col [],
3029 Int A []
3030 )
3031 {
3032 /* === Local variables ================================================== */
3033
3034 Int r ;
3035 Int c ;
3036 Int *rp ;
3037 Int *rp_end ;
3038 Int *cp ;
3039 Int *cp_end ;
3040
3041 /* === Dump the rows and columns of the matrix ========================== */
3042
3043 if (UMF_debug < 3)
3044 {
3045 return ;
3046 }
3047 DEBUG3 (("DUMP MATRIX:\n")) ;
3048 for (r = 0 ; r < n_row ; r++)
3049 {
3050 DEBUG3 (("Row "ID" alive? %d\n", r, ROW_IS_ALIVE (r))) ;
3051 if (ROW_IS_DEAD (r))
3052 {
3053 continue ;
3054 }
3055
3056 /* ------------------ */
3057 /* changed for UMFPACK: */
3058 DEBUG3 (("start "ID" length "ID" degree "ID" thickness "ID"\n",
3059 Row [r].start, Row [r].length, Row [r].shared1.degree,
3060 Row [r].thickness)) ;
3061 /* ------------------ */
3062
3063 rp = &A [Row [r].start] ;
3064 rp_end = rp + Row [r].length ;
3065 while (rp < rp_end)
3066 {
3067 c = *rp++ ;
3068 DEBUG4 ((" %d col "ID"\n", COL_IS_ALIVE (c), c)) ;
3069 }
3070 }
3071
3072 for (c = 0 ; c < n_col ; c++)
3073 {
3074 DEBUG3 (("Col "ID" alive? %d\n", c, COL_IS_ALIVE (c))) ;
3075 if (COL_IS_DEAD (c))
3076 {
3077 continue ;
3078 }
3079 /* ------------------ */
3080 /* changed for UMFPACK: */
3081 DEBUG3 (("start "ID" length "ID" shared1[thickness,parent] "ID
3082 " shared2 [order,score] "ID"\n", Col [c].start, Col [c].length,
3083 Col [c].shared1.thickness, Col [c].shared2.score));
3084 /* ------------------ */
3085 cp = &A [Col [c].start] ;
3086 cp_end = cp + Col [c].length ;
3087 while (cp < cp_end)
3088 {
3089 r = *cp++ ;
3090 DEBUG4 ((" %d row "ID"\n", ROW_IS_ALIVE (r), r)) ;
3091 }
3092
3093 /* ------------------ */
3094 /* added for UMFPACK: */
3095 DEBUG1 (("Col")) ;
3096 dump_super (c, Col, n_col) ;
3097 /* ------------------ */
3098
3099 }
3100 }
3101
3102 /* ------------------ */
3103 /* dump_super added for UMFPACK: */
dump_super(Int super_c,Colamd_Col Col[],Int n_col)3104 PRIVATE void dump_super
3105 (
3106 Int super_c,
3107 Colamd_Col Col [],
3108 Int n_col
3109 )
3110 {
3111 Int col, ncols ;
3112
3113 DEBUG1 ((" =[ ")) ;
3114 ncols = 0 ;
3115 for (col = super_c ; col != EMPTY ; col = Col [col].nextcol)
3116 {
3117 DEBUG1 ((" "ID, col)) ;
3118 ASSERT (col >= 0 && col < n_col) ;
3119 if (col != super_c)
3120 {
3121 ASSERT (COL_IS_DEAD (col)) ;
3122 }
3123 if (Col [col].nextcol == EMPTY)
3124 {
3125 ASSERT (col == Col [super_c].lastcol) ;
3126 }
3127 ncols++ ;
3128 ASSERT (ncols <= Col [super_c].shared1.thickness) ;
3129 }
3130 ASSERT (ncols == Col [super_c].shared1.thickness) ;
3131 DEBUG1 (("]\n")) ;
3132 }
3133 /* ------------------ */
3134
3135
3136 #endif /* NDEBUG */
3137