1 /* ========================================================================== */
2 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
3 /* ========================================================================== */
7 colamd: an approximate minimum degree column ordering algorithm,
8 for LU factorization of symmetric or unsymmetric matrices,
9 QR factorization, least squares, interior point methods for
10 linear programming problems, and other related problems.
12 symamd: an approximate minimum degree ordering algorithm for Cholesky
13 factorization of symmetric matrices.
17 Colamd computes a permutation Q such that the Cholesky factorization of
18 (AQ)'(AQ) has less fill-in and requires fewer floating point operations
19 than A'A. This also provides a good ordering for sparse partial
20 pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
21 factorization, and P is computed during numerical factorization via
22 conventional partial pivoting with row interchanges. Colamd is the
23 column ordering method used in SuperLU, part of the ScaLAPACK library.
24 It is also available as built-in function in MATLAB Version 6,
25 available from MathWorks, Inc. (http://www.mathworks.com). This
26 routine can be used in place of colmmd in MATLAB.
28 Symamd computes a permutation P of a symmetric matrix A such that the
29 Cholesky factorization of PAP' has less fill-in and requires fewer
30 floating point operations than A. Symamd constructs a matrix M such
31 that M'M has the same nonzero pattern of A, and then orders the columns
32 of M using colmmd. The column ordering of M is then returned as the
33 row and column ordering P of A.
37 The authors of the code itself are Stefan I. Larimore and Timothy A.
38 Davis (davis at cise.ufl.edu), University of Florida. The algorithm was
39 developed in collaboration with John Gilbert, Xerox PARC, and Esmond
40 Ng, Oak Ridge National Laboratory.
44 This work was supported by the National Science Foundation, under
45 grants DMS-9504974 and DMS-9803599.
47 Copyright and License:
49 Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
50 COLAMD is also available under alternate licenses, contact T. Davis
53 This library is free software; you can redistribute it and/or
54 modify it under the terms of the GNU Lesser General Public
55 License as published by the Free Software Foundation; either
56 version 2.1 of the License, or (at your option) any later version.
58 This library is distributed in the hope that it will be useful,
59 but WITHOUT ANY WARRANTY; without even the implied warranty of
60 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
61 Lesser General Public License for more details.
63 You should have received a copy of the GNU Lesser General Public
64 License along with this library; if not, write to the Free Software
65 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
68 Permission is hereby granted to use or copy this program under the
69 terms of the GNU LGPL, provided that the Copyright, this License,
70 and the Availability of the original version is retained on all copies.
71 User documentation of any code that uses this code or any modified
72 version of this code must cite the Copyright, this License, the
73 Availability note, and "Used by permission." Permission to modify
74 the code and to distribute modified code is granted, provided the
75 Copyright, this License, and the Availability note are retained,
76 and a notice that the code was modified is included.
80 The colamd/symamd library is available at
82 http://www.cise.ufl.edu/research/sparse/colamd/
84 This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
85 file. It requires the colamd.h file. It is required by the colamdmex.c
86 and symamdmex.c files, for the MATLAB interface to colamd and symamd.
87 Appears as ACM Algorithm 836.
89 See the ChangeLog file for changes since Version 1.0.
93 T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
94 minimum degree ordering algorithm, ACM Transactions on Mathematical
95 Software, vol. 30, no. 3., pp. 353-376, 2004.
97 T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
98 an approximate column minimum degree ordering algorithm, ACM
99 Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
104 /* ========================================================================== */
105 /* === Description of user-callable routines ================================ */
106 /* ========================================================================== */
108 /* COLAMD includes both int and UF_long versions of all its routines. The
109 * description below is for the int version. For UF_long, all int arguments
110 * become UF_long. UF_long is normally defined as long, except for WIN64.
112 ----------------------------------------------------------------------------
114 ----------------------------------------------------------------------------
119 size_t colamd_recommended (int nnz, int n_row, int n_col) ;
120 size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
125 Returns recommended value of Alen for use by colamd. Returns 0
126 if any input argument is negative. The use of this routine
127 is optional. Not needed for symamd, which dynamically allocates
130 Note that in v2.4 and earlier, these routines returned int or long.
131 They now return a value of type size_t.
133 Arguments (all input arguments):
135 int nnz ; Number of nonzeros in the matrix A. This must
136 be the same value as p [n_col] in the call to
137 colamd - otherwise you will get a wrong value
138 of the recommended memory to use.
140 int n_row ; Number of rows in the matrix A.
142 int n_col ; Number of columns in the matrix A.
144 ----------------------------------------------------------------------------
146 ----------------------------------------------------------------------------
151 colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
152 colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
156 Sets the default parameters. The use of this routine is optional.
160 double knobs [COLAMD_KNOBS] ; Output only.
162 NOTE: the meaning of the dense row/col knobs has changed in v2.4
164 knobs [0] and knobs [1] control dense row and col detection:
166 Colamd: rows with more than
167 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
168 entries are removed prior to ordering. Columns with more than
169 max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
170 entries are removed prior to
171 ordering, and placed last in the output column ordering.
173 Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
174 Rows and columns with more than
175 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
176 entries are removed prior to ordering, and placed last in the
179 COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
180 respectively, in colamd.h. Default values of these two knobs
181 are both 10. Currently, only knobs [0] and knobs [1] are
182 used, but future versions may use more knobs. If so, they will
183 be properly set to their defaults by the future version of
184 colamd_set_defaults, so that the code that calls colamd will
185 not need to change, assuming that you either use
186 colamd_set_defaults, or pass a (double *) NULL pointer as the
187 knobs array to colamd or symamd.
189 knobs [2]: aggressive absorption
191 knobs [COLAMD_AGGRESSIVE] controls whether or not to do
192 aggressive absorption during the ordering. Default is TRUE.
195 ----------------------------------------------------------------------------
197 ----------------------------------------------------------------------------
202 int colamd (int n_row, int n_col, int Alen, int *A, int *p,
203 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
204 UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
205 UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
206 UF_long stats [COLAMD_STATS]) ;
210 Computes a column ordering (Q) of A such that P(AQ)=LU or
211 (AQ)'AQ=LL' have less fill-in and require fewer floating point
212 operations than factorizing the unpermuted matrix A or A'A,
217 TRUE (1) if successful, FALSE (0) otherwise.
221 int n_row ; Input argument.
223 Number of rows in the matrix A.
224 Restriction: n_row >= 0.
225 Colamd returns FALSE if n_row is negative.
227 int n_col ; Input argument.
229 Number of columns in the matrix A.
230 Restriction: n_col >= 0.
231 Colamd returns FALSE if n_col is negative.
233 int Alen ; Input argument.
235 Restriction (see note):
236 Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
237 Colamd returns FALSE if these conditions are not met.
239 Note: this restriction makes an modest assumption regarding
240 the size of the two typedef's structures in colamd.h.
241 We do, however, guarantee that
243 Alen >= colamd_recommended (nnz, n_row, n_col)
245 will be sufficient. Note: the macro version does not check
246 for integer overflow, and thus is not recommended. Use
247 the colamd_recommended routine instead.
249 int A [Alen] ; Input argument, undefined on output.
251 A is an integer array of size Alen. Alen must be at least as
252 large as the bare minimum value given above, but this is very
253 low, and can result in excessive run time. For best
254 performance, we recommend that Alen be greater than or equal to
255 colamd_recommended (nnz, n_row, n_col), which adds
256 nnz/5 to the bare minimum value given above.
258 On input, the row indices of the entries in column c of the
259 matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices
260 in a given column c need not be in ascending order, and
261 duplicate row indices may be be present. However, colamd will
262 work a little faster if both of these conditions are met
263 (Colamd puts the matrix into this format, if it finds that the
264 the conditions are not met).
266 The matrix is 0-based. That is, rows are in the range 0 to
267 n_row-1, and columns are in the range 0 to n_col-1. Colamd
268 returns FALSE if any row index is out of range.
270 The contents of A are modified during ordering, and are
273 int p [n_col+1] ; Both input and output argument.
275 p is an integer array of size n_col+1. On input, it holds the
276 "pointers" for the column form of the matrix A. Column c of
277 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
278 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
279 for all c in the range 0 to n_col-1. The value p [n_col] is
280 thus the total number of entries in the pattern of the matrix A.
281 Colamd returns FALSE if these conditions are not met.
283 On output, if colamd returns TRUE, the array p holds the column
284 permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
285 the first column index in the new ordering, and p [n_col-1] is
286 the last. That is, p [k] = j means that column j of A is the
287 kth pivot column, in AQ, where k is in the range 0 to n_col-1
288 (p [0] = j means that column j of A is the first column in AQ).
290 If colamd returns FALSE, then no permutation is returned, and
291 p is undefined on output.
293 double knobs [COLAMD_KNOBS] ; Input argument.
295 See colamd_set_defaults for a description.
297 int stats [COLAMD_STATS] ; Output argument.
299 Statistics on the ordering, and error status.
300 See colamd.h for related definitions.
301 Colamd returns FALSE if stats is not present.
303 stats [0]: number of dense or empty rows ignored.
305 stats [1]: number of dense or empty columns ignored (and
306 ordered last in the output permutation p)
307 Note that a row can become "empty" if it
308 contains only "dense" and/or "empty" columns,
309 and similarly a column can become "empty" if it
310 only contains "dense" and/or "empty" rows.
312 stats [2]: number of garbage collections performed.
313 This can be excessively high if Alen is close
314 to the minimum required value.
316 stats [3]: status code. < 0 is an error code.
317 > 1 is a warning or notice.
319 0 OK. Each column of the input matrix contained
320 row indices in increasing order, with no
323 1 OK, but columns of input matrix were jumbled
324 (unsorted columns or duplicate entries). Colamd
325 had to do some extra work to sort the matrix
326 first and remove duplicate entries, but it
327 still was able to return a valid permutation
328 (return value of colamd was TRUE).
330 stats [4]: highest numbered column that
331 is unsorted or has duplicate
333 stats [5]: last seen duplicate or
335 stats [6]: number of duplicate or
336 unsorted row indices.
338 -1 A is a null pointer
340 -2 p is a null pointer
350 -5 number of nonzeros in matrix is negative
352 stats [4]: number of nonzeros, p [n_col]
360 stats [4]: required size
361 stats [5]: actual size (Alen)
363 -8 a column has a negative number of entries
365 stats [4]: column with < 0 entries
366 stats [5]: number of entries in col
368 -9 a row index is out of bounds
370 stats [4]: column with bad row index
371 stats [5]: bad row index
372 stats [6]: n_row, # of rows of matrx
374 -10 (unused; see symamd.c)
376 -999 (unused; see symamd.c)
378 Future versions may return more statistics in the stats array.
382 See http://www.cise.ufl.edu/research/sparse/colamd/example.c
383 for a complete example.
385 To order the columns of a 5-by-4 matrix with 11 nonzero entries in
386 the following nonzero pattern
394 with default knobs and no output statistics, do the following:
398 int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
399 int p [ ] = {0, 3, 5, 9, 11} ;
400 int stats [COLAMD_STATS] ;
401 colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
403 The permutation is returned in the array p, and A is destroyed.
405 ----------------------------------------------------------------------------
407 ----------------------------------------------------------------------------
412 int symamd (int n, int *A, int *p, int *perm,
413 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
414 void (*allocate) (size_t, size_t), void (*release) (void *)) ;
415 UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
416 double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
417 void (*allocate) (size_t, size_t), void (*release) (void *)) ;
421 The symamd routine computes an ordering P of a symmetric sparse
422 matrix A such that the Cholesky factorization PAP' = LL' remains
423 sparse. It is based on a column ordering of a matrix M constructed
424 so that the nonzero pattern of M'M is the same as A. The matrix A
425 is assumed to be symmetric; only the strictly lower triangular part
426 is accessed. You must pass your selected memory allocator (usually
427 calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
428 memory for the temporary matrix M.
432 TRUE (1) if successful, FALSE (0) otherwise.
436 int n ; Input argument.
438 Number of rows and columns in the symmetrix matrix A.
440 Symamd returns FALSE if n is negative.
442 int A [nnz] ; Input argument.
444 A is an integer array of size nnz, where nnz = p [n].
446 The row indices of the entries in column c of the matrix are
447 held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a
448 given column c need not be in ascending order, and duplicate
449 row indices may be present. However, symamd will run faster
450 if the columns are in sorted order with no duplicate entries.
452 The matrix is 0-based. That is, rows are in the range 0 to
453 n-1, and columns are in the range 0 to n-1. Symamd
454 returns FALSE if any row index is out of range.
456 The contents of A are not modified.
458 int p [n+1] ; Input argument.
460 p is an integer array of size n+1. On input, it holds the
461 "pointers" for the column form of the matrix A. Column c of
462 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
463 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
464 for all c in the range 0 to n-1. The value p [n] is
465 thus the total number of entries in the pattern of the matrix A.
466 Symamd returns FALSE if these conditions are not met.
468 The contents of p are not modified.
470 int perm [n+1] ; Output argument.
472 On output, if symamd returns TRUE, the array perm holds the
473 permutation P, where perm [0] is the first index in the new
474 ordering, and perm [n-1] is the last. That is, perm [k] = j
475 means that row and column j of A is the kth column in PAP',
476 where k is in the range 0 to n-1 (perm [0] = j means
477 that row and column j of A are the first row and column in
478 PAP'). The array is used as a workspace during the ordering,
479 which is why it must be of length n+1, not just n.
481 double knobs [COLAMD_KNOBS] ; Input argument.
483 See colamd_set_defaults for a description.
485 int stats [COLAMD_STATS] ; Output argument.
487 Statistics on the ordering, and error status.
488 See colamd.h for related definitions.
489 Symamd returns FALSE if stats is not present.
491 stats [0]: number of dense or empty row and columns ignored
492 (and ordered last in the output permutation
493 perm). Note that a row/column can become
494 "empty" if it contains only "dense" and/or
495 "empty" columns/rows.
497 stats [1]: (same as stats [0])
499 stats [2]: number of garbage collections performed.
501 stats [3]: status code. < 0 is an error code.
502 > 1 is a warning or notice.
504 0 OK. Each column of the input matrix contained
505 row indices in increasing order, with no
508 1 OK, but columns of input matrix were jumbled
509 (unsorted columns or duplicate entries). Symamd
510 had to do some extra work to sort the matrix
511 first and remove duplicate entries, but it
512 still was able to return a valid permutation
513 (return value of symamd was TRUE).
515 stats [4]: highest numbered column that
516 is unsorted or has duplicate
518 stats [5]: last seen duplicate or
520 stats [6]: number of duplicate or
521 unsorted row indices.
523 -1 A is a null pointer
525 -2 p is a null pointer
527 -3 (unused, see colamd.c)
533 -5 number of nonzeros in matrix is negative
535 stats [4]: # of nonzeros (p [n]).
543 -8 a column has a negative number of entries
545 stats [4]: column with < 0 entries
546 stats [5]: number of entries in col
548 -9 a row index is out of bounds
550 stats [4]: column with bad row index
551 stats [5]: bad row index
552 stats [6]: n_row, # of rows of matrx
554 -10 out of memory (unable to allocate temporary
555 workspace for M or count arrays using the
556 "allocate" routine passed into symamd).
558 Future versions may return more statistics in the stats array.
560 void * (*allocate) (size_t, size_t)
562 A pointer to a function providing memory allocation. The
563 allocated memory must be returned initialized to zero. For a
564 C application, this argument should normally be a pointer to
565 calloc. For a MATLAB mexFunction, the routine mxCalloc is
568 void (*release) (size_t, size_t)
570 A pointer to a function that frees memory allocated by the
571 memory allocation routine above. For a C application, this
572 argument should normally be a pointer to free. For a MATLAB
573 mexFunction, the routine mxFree is passed instead.
576 ----------------------------------------------------------------------------
578 ----------------------------------------------------------------------------
583 colamd_report (int stats [COLAMD_STATS]) ;
584 colamd_l_report (UF_long stats [COLAMD_STATS]) ;
588 Prints the error status and statistics recorded in the stats
589 array on the standard error output (for a standard C routine)
590 or on the MATLAB output (for a mexFunction).
594 int stats [COLAMD_STATS] ; Input only. Statistics from colamd.
597 ----------------------------------------------------------------------------
599 ----------------------------------------------------------------------------
604 symamd_report (int stats [COLAMD_STATS]) ;
605 symamd_l_report (UF_long stats [COLAMD_STATS]) ;
609 Prints the error status and statistics recorded in the stats
610 array on the standard error output (for a standard C routine)
611 or on the MATLAB output (for a mexFunction).
615 int stats [COLAMD_STATS] ; Input only. Statistics from symamd.
620 /* ========================================================================== */
621 /* === Scaffolding code definitions ======================================== */
622 /* ========================================================================== */
624 /* Ensure that debugging is turned off: */
629 /* turn on debugging by uncommenting the following line
634 Our "scaffolding code" philosophy: In our opinion, well-written library
635 code should keep its "debugging" code, and just normally have it turned off
636 by the compiler so as not to interfere with performance. This serves
639 (1) assertions act as comments to the reader, telling you what the code
640 expects at that point. All assertions will always be true (unless
641 there really is a bug, of course).
643 (2) leaving in the scaffolding code assists anyone who would like to modify
644 the code, or understand the algorithm (by reading the debugging output,
645 one can get a glimpse into what the code is doing).
647 (3) (gasp!) for actually finding bugs. This code has been heavily tested
648 and "should" be fully functional and bug-free ... but you never know...
650 The code will become outrageously slow when debugging is
651 enabled. To control the level of debugging output, set an environment
652 variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging,
653 you should see the following message on the standard output:
655 colamd: debug version, D = 1 (THIS WILL BE SLOW!)
657 or a similar message for symamd. If you don't, then debugging has not
662 /* ========================================================================== */
663 /* === Include files ======================================================== */
664 /* ========================================================================== */
672 #ifdef MATLAB_MEX_FILE
675 #endif /* MATLAB_MEX_FILE */
677 #if !defined (NPRINT) || !defined (NDEBUG)
682 #define NULL ((void *) 0)
686 /* ========================================================================== */
687 /* === int or UF_long ======================================================= */
688 /* ========================================================================== */
692 #include "UFconfig.h"
698 #define ID UF_long_id
699 #define Int_MAX UF_long_max
701 #define COLAMD_recommended colamd_l_recommended
702 #define COLAMD_set_defaults colamd_l_set_defaults
703 #define COLAMD_MAIN colamd_l
704 #define SYMAMD_MAIN symamd_l
705 #define COLAMD_report colamd_l_report
706 #define SYMAMD_report symamd_l_report
712 #define Int_MAX INT_MAX
714 #define COLAMD_recommended colamd_recommended
715 #define COLAMD_set_defaults colamd_set_defaults
716 #define COLAMD_MAIN colamd
717 #define SYMAMD_MAIN symamd
718 #define COLAMD_report colamd_report
719 #define SYMAMD_report symamd_report
723 /* ========================================================================== */
724 /* === Row and Column structures ============================================ */
725 /* ========================================================================== */
727 /* User code that makes use of the colamd/symamd routines need not directly */
728 /* reference these structures. They are used only for colamd_recommended. */
730 typedef struct Colamd_Col_struct
732 Int start ; /* index for A of first row in this column, or DEAD */
733 /* if column is dead */
734 Int length ; /* number of rows in this column */
737 Int thickness ; /* number of original columns represented by this */
738 /* col, if the column is alive */
739 Int parent ; /* parent in parent tree super-column structure, if */
740 /* the column is dead */
744 Int score ; /* the score used to maintain heap, if col is alive */
745 Int order ; /* pivot ordering of this column, if col is dead */
749 Int headhash ; /* head of a hash bucket, if col is at the head of */
751 Int hash ; /* hash value, if col is not in a degree list */
752 Int prev ; /* previous column in degree list, if col is in a */
753 /* degree list (but not at the head of a degree list) */
757 Int degree_next ; /* next column, if col is in a degree list */
758 Int hash_next ; /* next column, if col is in a hash list */
763 typedef struct Colamd_Row_struct
765 Int start ; /* index for A of first col in this row */
766 Int length ; /* number of principal columns in this row */
769 Int degree ; /* number of principal & non-principal columns in row */
770 Int p ; /* used as a row pointer in init_rows_cols () */
774 Int mark ; /* for computing set differences and marking dead rows*/
775 Int first_column ;/* first column in row (used in garbage collection) */
780 /* ========================================================================== */
781 /* === Definitions ========================================================== */
782 /* ========================================================================== */
784 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
786 #define PRIVATE static
788 #define DENSE_DEGREE(alpha,n) \
789 ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
791 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
792 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
794 #define ONES_COMPLEMENT(r) (-(r)-1)
796 /* -------------------------------------------------------------------------- */
797 /* Change for version 2.1: define TRUE and FALSE only if not yet defined */
798 /* -------------------------------------------------------------------------- */
808 /* -------------------------------------------------------------------------- */
812 /* Row and column status */
817 #define DEAD_PRINCIPAL (-1)
818 #define DEAD_NON_PRINCIPAL (-2)
820 /* Macros for row and column status update and checking. */
821 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
822 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
823 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
824 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
825 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
826 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
827 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
828 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
829 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
831 /* ========================================================================== */
832 /* === Colamd reporting mechanism =========================================== */
833 /* ========================================================================== */
835 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
836 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
837 #define INDEX(i) ((i)+1)
839 /* In C, matrices are 0-based and indices are reported as such in *_report */
843 /* All output goes through the PRINTF macro. */
844 #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
846 /* ========================================================================== */
847 /* === Prototypes of PRIVATE routines ======================================= */
848 /* ========================================================================== */
850 PRIVATE Int init_rows_cols
858 Int stats [COLAMD_STATS]
861 PRIVATE void init_scoring
869 double knobs [COLAMD_KNOBS],
875 PRIVATE Int find_ordering
890 PRIVATE void order_children
897 PRIVATE void detect_super_cols
912 PRIVATE Int garbage_collection
922 PRIVATE Int clear_mark
930 PRIVATE void print_report
933 Int stats [COLAMD_STATS]
936 /* ========================================================================== */
937 /* === Debugging prototypes and definitions ================================= */
938 /* ========================================================================== */
946 /* colamd_debug is the *ONLY* global variable, and is only */
947 /* present when debugging */
949 PRIVATE Int colamd_debug = 0 ; /* debug print level */
951 #define DEBUG0(params) { PRINTF (params) ; }
952 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
953 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
954 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
955 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
958 #ifdef MATLAB_MEX_FILE
959 #define ASSERT(expression) (mxAssert ((expression), ""))
961 #define ASSERT(expression) (assert (expression))
962 #endif /* MATLAB_MEX_FILE */
964 #define ASSERT xassert
967 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */
972 PRIVATE void debug_deg_lists
984 PRIVATE void debug_mark
992 PRIVATE void debug_matrix
1001 PRIVATE void debug_structures
1013 /* === No debugging ========================================================= */
1015 #define DEBUG0(params) ;
1016 #define DEBUG1(params) ;
1017 #define DEBUG2(params) ;
1018 #define DEBUG3(params) ;
1019 #define DEBUG4(params) ;
1021 #define ASSERT(expression)
1025 /* ========================================================================== */
1026 /* === USER-CALLABLE ROUTINES: ============================================== */
1027 /* ========================================================================== */
1029 /* ========================================================================== */
1030 /* === colamd_recommended =================================================== */
1031 /* ========================================================================== */
1034 The colamd_recommended routine returns the suggested size for Alen. This
1035 value has been determined to provide good balance between the number of
1036 garbage collections and the memory requirements for colamd. If any
1037 argument is negative, or if integer overflow occurs, a 0 is returned as an
1038 error condition. 2*nnz space is required for the row and column
1039 indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
1040 required for the Col and Row arrays, respectively, which are internal to
1041 colamd (roughly 6*n_col + 4*n_row). An additional n_col space is the
1042 minimal amount of "elbow room", and nnz/5 more space is recommended for
1043 run time efficiency.
1045 Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
1047 This function is not needed when using symamd.
1050 /* add two values of type size_t, and check for integer overflow */
1051 static size_t t_add (size_t a, size_t b, int *ok)
1053 (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
1054 return ((*ok) ? (a + b) : 0) ;
1057 /* compute a*k where k is a small integer, and check for integer overflow */
1058 static size_t t_mult (size_t a, size_t k, int *ok)
1061 for (i = 0 ; i < k ; i++)
1063 s = t_add (s, a, ok) ;
1068 /* size of the Col and Row structures */
1069 #define COLAMD_C(n_col,ok) \
1070 ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
1072 #define COLAMD_R(n_row,ok) \
1073 ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
1076 PUBLIC size_t COLAMD_recommended /* returns recommended value of Alen. */
1078 /* === Parameters ======================================================= */
1080 Int nnz, /* number of nonzeros in A */
1081 Int n_row, /* number of rows in A */
1082 Int n_col /* number of columns in A */
1087 if (nnz < 0 || n_row < 0 || n_col < 0)
1091 s = t_mult (nnz, 2, &ok) ; /* 2*nnz */
1092 c = COLAMD_C (n_col, &ok) ; /* size of column structures */
1093 r = COLAMD_R (n_row, &ok) ; /* size of row structures */
1094 s = t_add (s, c, &ok) ;
1095 s = t_add (s, r, &ok) ;
1096 s = t_add (s, n_col, &ok) ; /* elbow room */
1097 s = t_add (s, nnz/5, &ok) ; /* elbow room */
1098 ok = ok && (s < Int_MAX) ;
1099 return (ok ? s : 0) ;
1103 /* ========================================================================== */
1104 /* === colamd_set_defaults ================================================== */
1105 /* ========================================================================== */
1108 The colamd_set_defaults routine sets the default values of the user-
1109 controllable parameters for colamd and symamd:
1111 Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
1112 entries are removed prior to ordering. Columns with more than
1113 max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
1114 prior to ordering, and placed last in the output column ordering.
1116 Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
1117 entries are removed prior to ordering, and placed last in the
1120 knobs [0] dense row control
1122 knobs [1] dense column control
1124 knobs [2] if nonzero, do aggresive absorption
1126 knobs [3..19] unused, but future versions might use this
1130 PUBLIC void COLAMD_set_defaults
1132 /* === Parameters ======================================================= */
1134 double knobs [COLAMD_KNOBS] /* knob array */
1137 /* === Local variables ================================================== */
1143 return ; /* no knobs to initialize */
1145 for (i = 0 ; i < COLAMD_KNOBS ; i++)
1149 knobs [COLAMD_DENSE_ROW] = 10 ;
1150 knobs [COLAMD_DENSE_COL] = 10 ;
1151 knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
1155 /* ========================================================================== */
1156 /* === symamd =============================================================== */
1157 /* ========================================================================== */
1159 PUBLIC Int SYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */
1161 /* === Parameters ======================================================= */
1163 Int n, /* number of rows and columns of A */
1164 Int A [], /* row indices of A */
1165 Int p [], /* column pointers of A */
1166 Int perm [], /* output permutation, size n+1 */
1167 double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
1168 Int stats [COLAMD_STATS], /* output statistics and error codes */
1169 void * (*allocate) (size_t, size_t),
1170 /* pointer to calloc (ANSI C) or */
1171 /* mxCalloc (for MATLAB mexFunction) */
1172 void (*release) (void *)
1173 /* pointer to free (ANSI C) or */
1174 /* mxFree (for MATLAB mexFunction) */
1177 /* === Local variables ================================================== */
1179 Int *count ; /* length of each column of M, and col pointer*/
1180 Int *mark ; /* mark array for finding duplicate entries */
1181 Int *M ; /* row indices of matrix M */
1182 size_t Mlen ; /* length of M */
1183 Int n_row ; /* number of rows in M */
1184 Int nnz ; /* number of entries in A */
1185 Int i ; /* row index of A */
1186 Int j ; /* column index of A */
1187 Int k ; /* row index of M */
1188 Int mnz ; /* number of nonzeros in M */
1189 Int pp ; /* index into a column of A */
1190 Int last_row ; /* last row seen in the current column */
1191 Int length ; /* number of nonzeros in a column */
1193 double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */
1194 double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */
1197 colamd_get_debug ("symamd") ;
1200 /* === Check the input arguments ======================================== */
1204 DEBUG0 (("symamd: stats not present\n")) ;
1207 for (i = 0 ; i < COLAMD_STATS ; i++)
1211 stats [COLAMD_STATUS] = COLAMD_OK ;
1212 stats [COLAMD_INFO1] = -1 ;
1213 stats [COLAMD_INFO2] = -1 ;
1217 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1218 DEBUG0 (("symamd: A not present\n")) ;
1222 if (!p) /* p is not present */
1224 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1225 DEBUG0 (("symamd: p not present\n")) ;
1229 if (n < 0) /* n must be >= 0 */
1231 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1232 stats [COLAMD_INFO1] = n ;
1233 DEBUG0 (("symamd: n negative %d\n", n)) ;
1238 if (nnz < 0) /* nnz must be >= 0 */
1240 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1241 stats [COLAMD_INFO1] = nnz ;
1242 DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
1248 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1249 stats [COLAMD_INFO1] = p [0] ;
1250 DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
1254 /* === If no knobs, set default knobs =================================== */
1258 COLAMD_set_defaults (default_knobs) ;
1259 knobs = default_knobs ;
1262 /* === Allocate count and mark ========================================== */
1264 count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1267 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1268 DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
1272 mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1275 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1276 (*release) ((void *) count) ;
1277 DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
1281 /* === Compute column counts of M, check if A is valid ================== */
1283 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
1285 for (i = 0 ; i < n ; i++)
1290 for (j = 0 ; j < n ; j++)
1294 length = p [j+1] - p [j] ;
1297 /* column pointers must be non-decreasing */
1298 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1299 stats [COLAMD_INFO1] = j ;
1300 stats [COLAMD_INFO2] = length ;
1301 (*release) ((void *) count) ;
1302 (*release) ((void *) mark) ;
1303 DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
1307 for (pp = p [j] ; pp < p [j+1] ; pp++)
1310 if (i < 0 || i >= n)
1312 /* row index i, in column j, is out of bounds */
1313 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1314 stats [COLAMD_INFO1] = j ;
1315 stats [COLAMD_INFO2] = i ;
1316 stats [COLAMD_INFO3] = n ;
1317 (*release) ((void *) count) ;
1318 (*release) ((void *) mark) ;
1319 DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
1323 if (i <= last_row || mark [i] == j)
1325 /* row index is unsorted or repeated (or both), thus col */
1326 /* is jumbled. This is a notice, not an error condition. */
1327 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1328 stats [COLAMD_INFO1] = j ;
1329 stats [COLAMD_INFO2] = i ;
1330 (stats [COLAMD_INFO3]) ++ ;
1331 DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
1334 if (i > j && mark [i] != j)
1336 /* row k of M will contain column indices i and j */
1341 /* mark the row as having been seen in this column */
1348 /* v2.4: removed free(mark) */
1350 /* === Compute column pointers of M ===================================== */
1352 /* use output permutation, perm, for column pointers of M */
1354 for (j = 1 ; j <= n ; j++)
1356 perm [j] = perm [j-1] + count [j-1] ;
1358 for (j = 0 ; j < n ; j++)
1360 count [j] = perm [j] ;
1363 /* === Construct M ====================================================== */
1367 Mlen = COLAMD_recommended (mnz, n_row, n) ;
1368 M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
1369 DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
1370 n_row, n, mnz, (double) Mlen)) ;
1374 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1375 (*release) ((void *) count) ;
1376 (*release) ((void *) mark) ;
1377 DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
1383 if (stats [COLAMD_STATUS] == COLAMD_OK)
1386 for (j = 0 ; j < n ; j++)
1388 ASSERT (p [j+1] - p [j] >= 0) ;
1389 for (pp = p [j] ; pp < p [j+1] ; pp++)
1392 ASSERT (i >= 0 && i < n) ;
1395 /* row k of M contains column indices i and j */
1396 M [count [i]++] = k ;
1397 M [count [j]++] = k ;
1405 /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */
1406 DEBUG0 (("symamd: Duplicates in A.\n")) ;
1407 for (i = 0 ; i < n ; i++)
1411 for (j = 0 ; j < n ; j++)
1413 ASSERT (p [j+1] - p [j] >= 0) ;
1414 for (pp = p [j] ; pp < p [j+1] ; pp++)
1417 ASSERT (i >= 0 && i < n) ;
1418 if (i > j && mark [i] != j)
1420 /* row k of M contains column indices i and j */
1421 M [count [i]++] = k ;
1422 M [count [j]++] = k ;
1428 /* v2.4: free(mark) moved below */
1431 /* count and mark no longer needed */
1432 (*release) ((void *) count) ;
1433 (*release) ((void *) mark) ; /* v2.4: free (mark) moved here */
1434 ASSERT (k == n_row) ;
1436 /* === Adjust the knobs for M =========================================== */
1438 for (i = 0 ; i < COLAMD_KNOBS ; i++)
1440 cknobs [i] = knobs [i] ;
1443 /* there are no dense rows in M */
1444 cknobs [COLAMD_DENSE_ROW] = -1 ;
1445 cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
1447 /* === Order the columns of M =========================================== */
1449 /* v2.4: colamd cannot fail here, so the error check is removed */
1450 (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
1452 /* Note that the output permutation is now in perm */
1454 /* === get the statistics for symamd from colamd ======================== */
1456 /* a dense column in colamd means a dense row and col in symamd */
1457 stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
1459 /* === Free M =========================================================== */
1461 (*release) ((void *) M) ;
1462 DEBUG0 (("symamd: done.\n")) ;
1467 /* ========================================================================== */
1468 /* === colamd =============================================================== */
1469 /* ========================================================================== */
1472 The colamd routine computes a column ordering Q of a sparse matrix
1473 A such that the LU factorization P(AQ) = LU remains sparse, where P is
1474 selected via partial pivoting. The routine can also be viewed as
1475 providing a permutation Q such that the Cholesky factorization
1476 (AQ)'(AQ) = LL' remains sparse.
1479 PUBLIC Int COLAMD_MAIN /* returns TRUE if successful, FALSE otherwise*/
1481 /* === Parameters ======================================================= */
1483 Int n_row, /* number of rows in A */
1484 Int n_col, /* number of columns in A */
1485 Int Alen, /* length of A */
1486 Int A [], /* row indices of A */
1487 Int p [], /* pointers to columns in A */
1488 double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
1489 Int stats [COLAMD_STATS] /* output statistics and error codes */
1492 /* === Local variables ================================================== */
1494 Int i ; /* loop index */
1495 Int nnz ; /* nonzeros in A */
1496 size_t Row_size ; /* size of Row [], in integers */
1497 size_t Col_size ; /* size of Col [], in integers */
1498 size_t need ; /* minimum required length of A */
1499 Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */
1500 Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */
1501 Int n_col2 ; /* number of non-dense, non-empty columns */
1502 Int n_row2 ; /* number of non-dense, non-empty rows */
1503 Int ngarbage ; /* number of garbage collections performed */
1504 Int max_deg ; /* maximum row degree */
1505 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
1506 Int aggressive ; /* do aggressive absorption */
1510 colamd_get_debug ("colamd") ;
1513 /* === Check the input arguments ======================================== */
1517 DEBUG0 (("colamd: stats not present\n")) ;
1520 for (i = 0 ; i < COLAMD_STATS ; i++)
1524 stats [COLAMD_STATUS] = COLAMD_OK ;
1525 stats [COLAMD_INFO1] = -1 ;
1526 stats [COLAMD_INFO2] = -1 ;
1528 if (!A) /* A is not present */
1530 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1531 DEBUG0 (("colamd: A not present\n")) ;
1535 if (!p) /* p is not present */
1537 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1538 DEBUG0 (("colamd: p not present\n")) ;
1542 if (n_row < 0) /* n_row must be >= 0 */
1544 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
1545 stats [COLAMD_INFO1] = n_row ;
1546 DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
1550 if (n_col < 0) /* n_col must be >= 0 */
1552 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1553 stats [COLAMD_INFO1] = n_col ;
1554 DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
1559 if (nnz < 0) /* nnz must be >= 0 */
1561 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1562 stats [COLAMD_INFO1] = nnz ;
1563 DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
1569 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1570 stats [COLAMD_INFO1] = p [0] ;
1571 DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
1575 /* === If no knobs, set default knobs =================================== */
1579 COLAMD_set_defaults (default_knobs) ;
1580 knobs = default_knobs ;
1583 aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
1585 /* === Allocate the Row and Col arrays from array A ===================== */
1588 Col_size = COLAMD_C (n_col, &ok) ; /* size of Col array of structs */
1589 Row_size = COLAMD_R (n_row, &ok) ; /* size of Row array of structs */
1591 /* need = 2*nnz + n_col + Col_size + Row_size ; */
1592 need = t_mult (nnz, 2, &ok) ;
1593 need = t_add (need, n_col, &ok) ;
1594 need = t_add (need, Col_size, &ok) ;
1595 need = t_add (need, Row_size, &ok) ;
1597 if (!ok || need > (size_t) Alen || need > Int_MAX)
1599 /* not enough space in array A to perform the ordering */
1600 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
1601 stats [COLAMD_INFO1] = need ;
1602 stats [COLAMD_INFO2] = Alen ;
1603 DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
1607 Alen -= Col_size + Row_size ;
1608 Col = (Colamd_Col *) &A [Alen] ;
1609 Row = (Colamd_Row *) &A [Alen + Col_size] ;
1611 /* === Construct the row and column data structures ===================== */
1613 if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
1615 /* input matrix is invalid */
1616 DEBUG0 (("colamd: Matrix invalid\n")) ;
1620 /* === Initialize scores, kill dense rows/columns ======================= */
1622 init_scoring (n_row, n_col, Row, Col, A, p, knobs,
1623 &n_row2, &n_col2, &max_deg) ;
1625 /* === Order the supercolumns =========================================== */
1627 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
1628 n_col2, max_deg, 2*nnz, aggressive) ;
1630 /* === Order the non-principal columns ================================== */
1632 order_children (n_col, Col, p) ;
1634 /* === Return statistics in stats ======================================= */
1636 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
1637 stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
1638 stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
1639 DEBUG0 (("colamd: done.\n")) ;
1644 /* ========================================================================== */
1645 /* === colamd_report ======================================================== */
1646 /* ========================================================================== */
1648 PUBLIC void COLAMD_report
1650 Int stats [COLAMD_STATS]
1653 print_report ("colamd", stats) ;
1657 /* ========================================================================== */
1658 /* === symamd_report ======================================================== */
1659 /* ========================================================================== */
1661 PUBLIC void SYMAMD_report
1663 Int stats [COLAMD_STATS]
1666 print_report ("symamd", stats) ;
1671 /* ========================================================================== */
1672 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
1673 /* ========================================================================== */
1675 /* There are no user-callable routines beyond this point in the file */
1678 /* ========================================================================== */
1679 /* === init_rows_cols ======================================================= */
1680 /* ========================================================================== */
1683 Takes the column form of the matrix in A and creates the row form of the
1684 matrix. Also, row and column attributes are stored in the Col and Row
1685 structs. If the columns are un-sorted or contain duplicate row indices,
1686 this routine will also sort and remove duplicate row indices from the
1687 column form of the matrix. Returns FALSE if the matrix is invalid,
1688 TRUE otherwise. Not user-callable.
1691 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
1693 /* === Parameters ======================================================= */
1695 Int n_row, /* number of rows of A */
1696 Int n_col, /* number of columns of A */
1697 Colamd_Row Row [], /* of size n_row+1 */
1698 Colamd_Col Col [], /* of size n_col+1 */
1699 Int A [], /* row indices of A, of size Alen */
1700 Int p [], /* pointers to columns in A, of size n_col+1 */
1701 Int stats [COLAMD_STATS] /* colamd statistics */
1704 /* === Local variables ================================================== */
1706 Int col ; /* a column index */
1707 Int row ; /* a row index */
1708 Int *cp ; /* a column pointer */
1709 Int *cp_end ; /* a pointer to the end of a column */
1710 Int *rp ; /* a row pointer */
1711 Int *rp_end ; /* a pointer to the end of a row */
1712 Int last_row ; /* previous row */
1714 /* === Initialize columns, and check column pointers ==================== */
1716 for (col = 0 ; col < n_col ; col++)
1718 Col [col].start = p [col] ;
1719 Col [col].length = p [col+1] - p [col] ;
1721 if (Col [col].length < 0)
1723 /* column pointers must be non-decreasing */
1724 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1725 stats [COLAMD_INFO1] = col ;
1726 stats [COLAMD_INFO2] = Col [col].length ;
1727 DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
1731 Col [col].shared1.thickness = 1 ;
1732 Col [col].shared2.score = 0 ;
1733 Col [col].shared3.prev = EMPTY ;
1734 Col [col].shared4.degree_next = EMPTY ;
1737 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
1739 /* === Scan columns, compute row degrees, and check row indices ========= */
1741 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
1743 for (row = 0 ; row < n_row ; row++)
1745 Row [row].length = 0 ;
1746 Row [row].shared2.mark = -1 ;
1749 for (col = 0 ; col < n_col ; col++)
1754 cp_end = &A [p [col+1]] ;
1760 /* make sure row indices within range */
1761 if (row < 0 || row >= n_row)
1763 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1764 stats [COLAMD_INFO1] = col ;
1765 stats [COLAMD_INFO2] = row ;
1766 stats [COLAMD_INFO3] = n_row ;
1767 DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
1771 if (row <= last_row || Row [row].shared2.mark == col)
1773 /* row index are unsorted or repeated (or both), thus col */
1774 /* is jumbled. This is a notice, not an error condition. */
1775 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1776 stats [COLAMD_INFO1] = col ;
1777 stats [COLAMD_INFO2] = row ;
1778 (stats [COLAMD_INFO3]) ++ ;
1779 DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
1782 if (Row [row].shared2.mark != col)
1784 Row [row].length++ ;
1788 /* this is a repeated entry in the column, */
1789 /* it will be removed */
1790 Col [col].length-- ;
1793 /* mark the row as having been seen in this column */
1794 Row [row].shared2.mark = col ;
1800 /* === Compute row pointers ============================================= */
1802 /* row form of the matrix starts directly after the column */
1803 /* form of matrix in A */
1804 Row [0].start = p [n_col] ;
1805 Row [0].shared1.p = Row [0].start ;
1806 Row [0].shared2.mark = -1 ;
1807 for (row = 1 ; row < n_row ; row++)
1809 Row [row].start = Row [row-1].start + Row [row-1].length ;
1810 Row [row].shared1.p = Row [row].start ;
1811 Row [row].shared2.mark = -1 ;
1814 /* === Create row form ================================================== */
1816 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1818 /* if cols jumbled, watch for repeated row indices */
1819 for (col = 0 ; col < n_col ; col++)
1822 cp_end = &A [p [col+1]] ;
1826 if (Row [row].shared2.mark != col)
1828 A [(Row [row].shared1.p)++] = col ;
1829 Row [row].shared2.mark = col ;
1836 /* if cols not jumbled, we don't need the mark (this is faster) */
1837 for (col = 0 ; col < n_col ; col++)
1840 cp_end = &A [p [col+1]] ;
1843 A [(Row [*cp++].shared1.p)++] = col ;
1848 /* === Clear the row marks and set row degrees ========================== */
1850 for (row = 0 ; row < n_row ; row++)
1852 Row [row].shared2.mark = 0 ;
1853 Row [row].shared1.degree = Row [row].length ;
1856 /* === See if we need to re-create columns ============================== */
1858 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1860 DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
1863 /* make sure column lengths are correct */
1864 for (col = 0 ; col < n_col ; col++)
1866 p [col] = Col [col].length ;
1868 for (row = 0 ; row < n_row ; row++)
1870 rp = &A [Row [row].start] ;
1871 rp_end = rp + Row [row].length ;
1877 for (col = 0 ; col < n_col ; col++)
1879 ASSERT (p [col] == 0) ;
1881 /* now p is all zero (different than when debugging is turned off) */
1884 /* === Compute col pointers ========================================= */
1886 /* col form of the matrix starts at A [0]. */
1887 /* Note, we may have a gap between the col form and the row */
1888 /* form if there were duplicate entries, if so, it will be */
1889 /* removed upon the first garbage collection */
1891 p [0] = Col [0].start ;
1892 for (col = 1 ; col < n_col ; col++)
1894 /* note that the lengths here are for pruned columns, i.e. */
1895 /* no duplicate row indices will exist for these columns */
1896 Col [col].start = Col [col-1].start + Col [col-1].length ;
1897 p [col] = Col [col].start ;
1900 /* === Re-create col form =========================================== */
1902 for (row = 0 ; row < n_row ; row++)
1904 rp = &A [Row [row].start] ;
1905 rp_end = rp + Row [row].length ;
1908 A [(p [*rp++])++] = row ;
1913 /* === Done. Matrix is not (or no longer) jumbled ====================== */
1919 /* ========================================================================== */
1920 /* === init_scoring ========================================================= */
1921 /* ========================================================================== */
1924 Kills dense or empty columns and rows, calculates an initial score for
1925 each column, and places all columns in the degree lists. Not user-callable.
1928 PRIVATE void init_scoring
1930 /* === Parameters ======================================================= */
1932 Int n_row, /* number of rows of A */
1933 Int n_col, /* number of columns of A */
1934 Colamd_Row Row [], /* of size n_row+1 */
1935 Colamd_Col Col [], /* of size n_col+1 */
1936 Int A [], /* column form and row form of A */
1937 Int head [], /* of size n_col+1 */
1938 double knobs [COLAMD_KNOBS],/* parameters */
1939 Int *p_n_row2, /* number of non-dense, non-empty rows */
1940 Int *p_n_col2, /* number of non-dense, non-empty columns */
1941 Int *p_max_deg /* maximum row degree */
1944 /* === Local variables ================================================== */
1946 Int c ; /* a column index */
1947 Int r, row ; /* a row index */
1948 Int *cp ; /* a column pointer */
1949 Int deg ; /* degree of a row or column */
1950 Int *cp_end ; /* a pointer to the end of a column */
1951 Int *new_cp ; /* new column pointer */
1952 Int col_length ; /* length of pruned column */
1953 Int score ; /* current column score */
1954 Int n_col2 ; /* number of non-dense, non-empty columns */
1955 Int n_row2 ; /* number of non-dense, non-empty rows */
1956 Int dense_row_count ; /* remove rows with more entries than this */
1957 Int dense_col_count ; /* remove cols with more entries than this */
1958 Int min_score ; /* smallest column score */
1959 Int max_deg ; /* maximum row degree */
1960 Int next_col ; /* Used to add to degree list.*/
1963 Int debug_count ; /* debug only. */
1966 /* === Extract knobs ==================================================== */
1968 /* Note: if knobs contains a NaN, this is undefined: */
1969 if (knobs [COLAMD_DENSE_ROW] < 0)
1971 /* only remove completely dense rows */
1972 dense_row_count = n_col-1 ;
1976 dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
1978 if (knobs [COLAMD_DENSE_COL] < 0)
1980 /* only remove completely dense columns */
1981 dense_col_count = n_row-1 ;
1986 DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
1989 DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
1994 /* === Kill empty columns =============================================== */
1996 /* Put the empty columns at the end in their natural order, so that LU */
1997 /* factorization can proceed as far as possible. */
1998 for (c = n_col-1 ; c >= 0 ; c--)
2000 deg = Col [c].length ;
2003 /* this is a empty column, kill and order it last */
2004 Col [c].shared2.order = --n_col2 ;
2005 KILL_PRINCIPAL_COL (c) ;
2008 DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
2010 /* === Kill dense columns =============================================== */
2012 /* Put the dense columns at the end, in their natural order */
2013 for (c = n_col-1 ; c >= 0 ; c--)
2015 /* skip any dead columns */
2016 if (COL_IS_DEAD (c))
2020 deg = Col [c].length ;
2021 if (deg > dense_col_count)
2023 /* this is a dense column, kill and order it last */
2024 Col [c].shared2.order = --n_col2 ;
2025 /* decrement the row degrees */
2026 cp = &A [Col [c].start] ;
2027 cp_end = cp + Col [c].length ;
2030 Row [*cp++].shared1.degree-- ;
2032 KILL_PRINCIPAL_COL (c) ;
2035 DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
2037 /* === Kill dense and empty rows ======================================== */
2039 for (r = 0 ; r < n_row ; r++)
2041 deg = Row [r].shared1.degree ;
2042 ASSERT (deg >= 0 && deg <= n_col) ;
2043 if (deg > dense_row_count || deg == 0)
2045 /* kill a dense or empty row */
2051 /* keep track of max degree of remaining rows */
2052 max_deg = MAX (max_deg, deg) ;
2055 DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
2057 /* === Compute initial column scores ==================================== */
2059 /* At this point the row degrees are accurate. They reflect the number */
2060 /* of "live" (non-dense) columns in each row. No empty rows exist. */
2061 /* Some "live" columns may contain only dead rows, however. These are */
2062 /* pruned in the code below. */
2064 /* now find the initial matlab score for each column */
2065 for (c = n_col-1 ; c >= 0 ; c--)
2067 /* skip dead column */
2068 if (COL_IS_DEAD (c))
2073 cp = &A [Col [c].start] ;
2075 cp_end = cp + Col [c].length ;
2081 if (ROW_IS_DEAD (row))
2085 /* compact the column */
2087 /* add row's external degree */
2088 score += Row [row].shared1.degree - 1 ;
2089 /* guard against integer overflow */
2090 score = MIN (score, n_col) ;
2092 /* determine pruned column length */
2093 col_length = (Int) (new_cp - &A [Col [c].start]) ;
2094 if (col_length == 0)
2096 /* a newly-made null column (all rows in this col are "dense" */
2097 /* and have already been killed) */
2098 DEBUG2 (("Newly null killed: %d\n", c)) ;
2099 Col [c].shared2.order = --n_col2 ;
2100 KILL_PRINCIPAL_COL (c) ;
2104 /* set column length and set score */
2105 ASSERT (score >= 0) ;
2106 ASSERT (score <= n_col) ;
2107 Col [c].length = col_length ;
2108 Col [c].shared2.score = score ;
2111 DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
2114 /* At this point, all empty rows and columns are dead. All live columns */
2115 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
2116 /* yet). Rows may contain dead columns, but all live rows contain at */
2117 /* least one live column. */
2120 debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
2123 /* === Initialize degree lists ========================================== */
2129 /* clear the hash buckets */
2130 for (c = 0 ; c <= n_col ; c++)
2135 /* place in reverse order, so low column indices are at the front */
2136 /* of the lists. This is to encourage natural tie-breaking */
2137 for (c = n_col-1 ; c >= 0 ; c--)
2139 /* only add principal columns to degree lists */
2140 if (COL_IS_ALIVE (c))
2142 DEBUG4 (("place %d score %d minscore %d ncol %d\n",
2143 c, Col [c].shared2.score, min_score, n_col)) ;
2145 /* === Add columns score to DList =============================== */
2147 score = Col [c].shared2.score ;
2149 ASSERT (min_score >= 0) ;
2150 ASSERT (min_score <= n_col) ;
2151 ASSERT (score >= 0) ;
2152 ASSERT (score <= n_col) ;
2153 ASSERT (head [score] >= EMPTY) ;
2155 /* now add this column to dList at proper score location */
2156 next_col = head [score] ;
2157 Col [c].shared3.prev = EMPTY ;
2158 Col [c].shared4.degree_next = next_col ;
2160 /* if there already was a column with the same score, set its */
2161 /* previous pointer to this new column */
2162 if (next_col != EMPTY)
2164 Col [next_col].shared3.prev = c ;
2168 /* see if this score is less than current min */
2169 min_score = MIN (min_score, score) ;
2179 DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
2180 debug_count, n_col, n_col-debug_count)) ;
2181 ASSERT (debug_count == n_col2) ;
2182 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
2185 /* === Return number of remaining columns, and max row degree =========== */
2187 *p_n_col2 = n_col2 ;
2188 *p_n_row2 = n_row2 ;
2189 *p_max_deg = max_deg ;
2193 /* ========================================================================== */
2194 /* === find_ordering ======================================================== */
2195 /* ========================================================================== */
2198 Order the principal columns of the supercolumn form of the matrix
2199 (no supercolumns on input). Uses a minimum approximate column minimum
2200 degree ordering method. Not user-callable.
2203 PRIVATE Int find_ordering /* return the number of garbage collections */
2205 /* === Parameters ======================================================= */
2207 Int n_row, /* number of rows of A */
2208 Int n_col, /* number of columns of A */
2209 Int Alen, /* size of A, 2*nnz + n_col or larger */
2210 Colamd_Row Row [], /* of size n_row+1 */
2211 Colamd_Col Col [], /* of size n_col+1 */
2212 Int A [], /* column form and row form of A */
2213 Int head [], /* of size n_col+1 */
2214 Int n_col2, /* Remaining columns to order */
2215 Int max_deg, /* Maximum row degree */
2216 Int pfree, /* index of first free slot (2*nnz on entry) */
2220 /* === Local variables ================================================== */
2222 Int k ; /* current pivot ordering step */
2223 Int pivot_col ; /* current pivot column */
2224 Int *cp ; /* a column pointer */
2225 Int *rp ; /* a row pointer */
2226 Int pivot_row ; /* current pivot row */
2227 Int *new_cp ; /* modified column pointer */
2228 Int *new_rp ; /* modified row pointer */
2229 Int pivot_row_start ; /* pointer to start of pivot row */
2230 Int pivot_row_degree ; /* number of columns in pivot row */
2231 Int pivot_row_length ; /* number of supercolumns in pivot row */
2232 Int pivot_col_score ; /* score of pivot column */
2233 Int needed_memory ; /* free space needed for pivot row */
2234 Int *cp_end ; /* pointer to the end of a column */
2235 Int *rp_end ; /* pointer to the end of a row */
2236 Int row ; /* a row index */
2237 Int col ; /* a column index */
2238 Int max_score ; /* maximum possible score */
2239 Int cur_score ; /* score of current column */
2240 unsigned Int hash ; /* hash value for supernode detection */
2241 Int head_column ; /* head of hash bucket */
2242 Int first_col ; /* first column in hash bucket */
2243 Int tag_mark ; /* marker value for mark array */
2244 Int row_mark ; /* Row [row].shared2.mark */
2245 Int set_difference ; /* set difference size of row with pivot row */
2246 Int min_score ; /* smallest column score */
2247 Int col_thickness ; /* "thickness" (no. of columns in a supercol) */
2248 Int max_mark ; /* maximum value of tag_mark */
2249 Int pivot_col_thickness ; /* number of columns represented by pivot col */
2250 Int prev_col ; /* Used by Dlist operations. */
2251 Int next_col ; /* Used by Dlist operations. */
2252 Int ngarbage ; /* number of garbage collections performed */
2255 Int debug_d ; /* debug loop counter */
2256 Int debug_step = 0 ; /* debug loop counter */
2259 /* === Initialization and clear mark ==================================== */
2261 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */
2262 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2265 DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
2267 /* === Order the columns ================================================ */
2269 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
2273 if (debug_step % 100 == 0)
2275 DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2279 DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2282 debug_deg_lists (n_row, n_col, Row, Col, head,
2283 min_score, n_col2-k, max_deg) ;
2284 debug_matrix (n_row, n_col, Row, Col, A) ;
2287 /* === Select pivot column, and order it ============================ */
2289 /* make sure degree list isn't empty */
2290 ASSERT (min_score >= 0) ;
2291 ASSERT (min_score <= n_col) ;
2292 ASSERT (head [min_score] >= EMPTY) ;
2295 for (debug_d = 0 ; debug_d < min_score ; debug_d++)
2297 ASSERT (head [debug_d] == EMPTY) ;
2301 /* get pivot column from head of minimum degree list */
2302 while (head [min_score] == EMPTY && min_score < n_col)
2306 pivot_col = head [min_score] ;
2307 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
2308 next_col = Col [pivot_col].shared4.degree_next ;
2309 head [min_score] = next_col ;
2310 if (next_col != EMPTY)
2312 Col [next_col].shared3.prev = EMPTY ;
2315 ASSERT (COL_IS_ALIVE (pivot_col)) ;
2317 /* remember score for defrag check */
2318 pivot_col_score = Col [pivot_col].shared2.score ;
2320 /* the pivot column is the kth column in the pivot order */
2321 Col [pivot_col].shared2.order = k ;
2323 /* increment order count by column thickness */
2324 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
2325 k += pivot_col_thickness ;
2326 ASSERT (pivot_col_thickness > 0) ;
2327 DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
2329 /* === Garbage_collection, if necessary ============================= */
2331 needed_memory = MIN (pivot_col_score, n_col - k) ;
2332 if (pfree + needed_memory >= Alen)
2334 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
2336 /* after garbage collection we will have enough */
2337 ASSERT (pfree + needed_memory < Alen) ;
2338 /* garbage collection has wiped out the Row[].shared2.mark array */
2339 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2342 debug_matrix (n_row, n_col, Row, Col, A) ;
2346 /* === Compute pivot row pattern ==================================== */
2348 /* get starting location for this new merged row */
2349 pivot_row_start = pfree ;
2351 /* initialize new row counts to zero */
2352 pivot_row_degree = 0 ;
2354 /* tag pivot column as having been visited so it isn't included */
2355 /* in merged pivot row */
2356 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
2358 /* pivot row is the union of all rows in the pivot column pattern */
2359 cp = &A [Col [pivot_col].start] ;
2360 cp_end = cp + Col [pivot_col].length ;
2365 DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
2366 /* skip if row is dead */
2367 if (ROW_IS_ALIVE (row))
2369 rp = &A [Row [row].start] ;
2370 rp_end = rp + Row [row].length ;
2375 /* add the column, if alive and untagged */
2376 col_thickness = Col [col].shared1.thickness ;
2377 if (col_thickness > 0 && COL_IS_ALIVE (col))
2379 /* tag column in pivot row */
2380 Col [col].shared1.thickness = -col_thickness ;
2381 ASSERT (pfree < Alen) ;
2382 /* place column in pivot row */
2384 pivot_row_degree += col_thickness ;
2390 /* clear tag on pivot column */
2391 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
2392 max_deg = MAX (max_deg, pivot_row_degree) ;
2395 DEBUG3 (("check2\n")) ;
2396 debug_mark (n_row, Row, tag_mark, max_mark) ;
2399 /* === Kill all rows used to construct pivot row ==================== */
2401 /* also kill pivot row, temporarily */
2402 cp = &A [Col [pivot_col].start] ;
2403 cp_end = cp + Col [pivot_col].length ;
2406 /* may be killing an already dead row */
2408 DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
2412 /* === Select a row index to use as the new pivot row =============== */
2414 pivot_row_length = pfree - pivot_row_start ;
2415 if (pivot_row_length > 0)
2417 /* pick the "pivot" row arbitrarily (first row in col) */
2418 pivot_row = A [Col [pivot_col].start] ;
2419 DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
2423 /* there is no pivot row, since it is of zero length */
2425 ASSERT (pivot_row_length == 0) ;
2427 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
2429 /* === Approximate degree computation =============================== */
2431 /* Here begins the computation of the approximate degree. The column */
2432 /* score is the sum of the pivot row "length", plus the size of the */
2433 /* set differences of each row in the column minus the pattern of the */
2434 /* pivot row itself. The column ("thickness") itself is also */
2435 /* excluded from the column score (we thus use an approximate */
2436 /* external degree). */
2438 /* The time taken by the following code (compute set differences, and */
2439 /* add them up) is proportional to the size of the data structure */
2440 /* being scanned - that is, the sum of the sizes of each column in */
2441 /* the pivot row. Thus, the amortized time to compute a column score */
2442 /* is proportional to the size of that column (where size, in this */
2443 /* context, is the column "length", or the number of row indices */
2444 /* in that column). The number of row indices in a column is */
2445 /* monotonically non-decreasing, from the length of the original */
2446 /* column on input to colamd. */
2448 /* === Compute set differences ====================================== */
2450 DEBUG3 (("** Computing set differences phase. **\n")) ;
2452 /* pivot row is currently dead - it will be revived later. */
2454 DEBUG3 (("Pivot row: ")) ;
2455 /* for each column in pivot row */
2456 rp = &A [pivot_row_start] ;
2457 rp_end = rp + pivot_row_length ;
2461 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2462 DEBUG3 (("Col: %d\n", col)) ;
2464 /* clear tags used to construct pivot row pattern */
2465 col_thickness = -Col [col].shared1.thickness ;
2466 ASSERT (col_thickness > 0) ;
2467 Col [col].shared1.thickness = col_thickness ;
2469 /* === Remove column from degree list =========================== */
2471 cur_score = Col [col].shared2.score ;
2472 prev_col = Col [col].shared3.prev ;
2473 next_col = Col [col].shared4.degree_next ;
2474 ASSERT (cur_score >= 0) ;
2475 ASSERT (cur_score <= n_col) ;
2476 ASSERT (cur_score >= EMPTY) ;
2477 if (prev_col == EMPTY)
2479 head [cur_score] = next_col ;
2483 Col [prev_col].shared4.degree_next = next_col ;
2485 if (next_col != EMPTY)
2487 Col [next_col].shared3.prev = prev_col ;
2490 /* === Scan the column ========================================== */
2492 cp = &A [Col [col].start] ;
2493 cp_end = cp + Col [col].length ;
2498 row_mark = Row [row].shared2.mark ;
2500 if (ROW_IS_MARKED_DEAD (row_mark))
2504 ASSERT (row != pivot_row) ;
2505 set_difference = row_mark - tag_mark ;
2506 /* check if the row has been seen yet */
2507 if (set_difference < 0)
2509 ASSERT (Row [row].shared1.degree <= max_deg) ;
2510 set_difference = Row [row].shared1.degree ;
2512 /* subtract column thickness from this row's set difference */
2513 set_difference -= col_thickness ;
2514 ASSERT (set_difference >= 0) ;
2515 /* absorb this row if the set difference becomes zero */
2516 if (set_difference == 0 && aggressive)
2518 DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
2523 /* save the new mark */
2524 Row [row].shared2.mark = set_difference + tag_mark ;
2530 debug_deg_lists (n_row, n_col, Row, Col, head,
2531 min_score, n_col2-k-pivot_row_degree, max_deg) ;
2534 /* === Add up set differences for each column ======================= */
2536 DEBUG3 (("** Adding set differences phase. **\n")) ;
2538 /* for each column in pivot row */
2539 rp = &A [pivot_row_start] ;
2540 rp_end = rp + pivot_row_length ;
2545 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2548 cp = &A [Col [col].start] ;
2549 /* compact the column */
2551 cp_end = cp + Col [col].length ;
2553 DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
2559 ASSERT(row >= 0 && row < n_row) ;
2560 row_mark = Row [row].shared2.mark ;
2562 if (ROW_IS_MARKED_DEAD (row_mark))
2564 DEBUG4 ((" Row %d, dead\n", row)) ;
2567 DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
2568 ASSERT (row_mark >= tag_mark) ;
2569 /* compact the column */
2571 /* compute hash function */
2573 /* add set difference */
2574 cur_score += row_mark - tag_mark ;
2575 /* integer overflow... */
2576 cur_score = MIN (cur_score, n_col) ;
2579 /* recompute the column's length */
2580 Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
2582 /* === Further mass elimination ================================= */
2584 if (Col [col].length == 0)
2586 DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
2587 /* nothing left but the pivot row in this column */
2588 KILL_PRINCIPAL_COL (col) ;
2589 pivot_row_degree -= Col [col].shared1.thickness ;
2590 ASSERT (pivot_row_degree >= 0) ;
2592 Col [col].shared2.order = k ;
2593 /* increment order count by column thickness */
2594 k += Col [col].shared1.thickness ;
2598 /* === Prepare for supercolumn detection ==================== */
2600 DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
2602 /* save score so far */
2603 Col [col].shared2.score = cur_score ;
2605 /* add column to hash table, for supercolumn detection */
2608 DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
2609 ASSERT (((Int) hash) <= n_col) ;
2611 head_column = head [hash] ;
2612 if (head_column > EMPTY)
2614 /* degree list "hash" is non-empty, use prev (shared3) of */
2615 /* first column in degree list as head of hash bucket */
2616 first_col = Col [head_column].shared3.headhash ;
2617 Col [head_column].shared3.headhash = col ;
2621 /* degree list "hash" is empty, use head as hash bucket */
2622 first_col = - (head_column + 2) ;
2623 head [hash] = - (col + 2) ;
2625 Col [col].shared4.hash_next = first_col ;
2627 /* save hash function in Col [col].shared3.hash */
2628 Col [col].shared3.hash = (Int) hash ;
2629 ASSERT (COL_IS_ALIVE (col)) ;
2633 /* The approximate external column degree is now computed. */
2635 /* === Supercolumn detection ======================================== */
2637 DEBUG3 (("** Supercolumn detection phase. **\n")) ;
2645 Col, A, head, pivot_row_start, pivot_row_length) ;
2647 /* === Kill the pivotal column ====================================== */
2649 KILL_PRINCIPAL_COL (pivot_col) ;
2651 /* === Clear mark =================================================== */
2653 tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
2656 DEBUG3 (("check3\n")) ;
2657 debug_mark (n_row, Row, tag_mark, max_mark) ;
2660 /* === Finalize the new pivot row, and column scores ================ */
2662 DEBUG3 (("** Finalize scores phase. **\n")) ;
2664 /* for each column in pivot row */
2665 rp = &A [pivot_row_start] ;
2666 /* compact the pivot row */
2668 rp_end = rp + pivot_row_length ;
2672 /* skip dead columns */
2673 if (COL_IS_DEAD (col))
2678 /* add new pivot row to column */
2679 A [Col [col].start + (Col [col].length++)] = pivot_row ;
2681 /* retrieve score so far and add on pivot row's degree. */
2682 /* (we wait until here for this in case the pivot */
2683 /* row's degree was reduced due to mass elimination). */
2684 cur_score = Col [col].shared2.score + pivot_row_degree ;
2686 /* calculate the max possible score as the number of */
2687 /* external columns minus the 'k' value minus the */
2688 /* columns thickness */
2689 max_score = n_col - k - Col [col].shared1.thickness ;
2691 /* make the score the external degree of the union-of-rows */
2692 cur_score -= Col [col].shared1.thickness ;
2694 /* make sure score is less or equal than the max score */
2695 cur_score = MIN (cur_score, max_score) ;
2696 ASSERT (cur_score >= 0) ;
2698 /* store updated score */
2699 Col [col].shared2.score = cur_score ;
2701 /* === Place column back in degree list ========================= */
2703 ASSERT (min_score >= 0) ;
2704 ASSERT (min_score <= n_col) ;
2705 ASSERT (cur_score >= 0) ;
2706 ASSERT (cur_score <= n_col) ;
2707 ASSERT (head [cur_score] >= EMPTY) ;
2708 next_col = head [cur_score] ;
2709 Col [col].shared4.degree_next = next_col ;
2710 Col [col].shared3.prev = EMPTY ;
2711 if (next_col != EMPTY)
2713 Col [next_col].shared3.prev = col ;
2715 head [cur_score] = col ;
2717 /* see if this score is less than current min */
2718 min_score = MIN (min_score, cur_score) ;
2723 debug_deg_lists (n_row, n_col, Row, Col, head,
2724 min_score, n_col2-k, max_deg) ;
2727 /* === Resurrect the new pivot row ================================== */
2729 if (pivot_row_degree > 0)
2731 /* update pivot row length to reflect any cols that were killed */
2732 /* during super-col detection and mass elimination */
2733 Row [pivot_row].start = pivot_row_start ;
2734 Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
2735 ASSERT (Row [pivot_row].length > 0) ;
2736 Row [pivot_row].shared1.degree = pivot_row_degree ;
2737 Row [pivot_row].shared2.mark = 0 ;
2738 /* pivot row is no longer dead */
2740 DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
2741 pivot_row, pivot_row_degree)) ;
2745 /* === All principal columns have now been ordered ====================== */
2751 /* ========================================================================== */
2752 /* === order_children ======================================================= */
2753 /* ========================================================================== */
2756 The find_ordering routine has ordered all of the principal columns (the
2757 representatives of the supercolumns). The non-principal columns have not
2758 yet been ordered. This routine orders those columns by walking up the
2759 parent tree (a column is a child of the column which absorbed it). The
2760 final permutation vector is then placed in p [0 ... n_col-1], with p [0]
2761 being the first column, and p [n_col-1] being the last. It doesn't look
2762 like it at first glance, but be assured that this routine takes time linear
2763 in the number of columns. Although not immediately obvious, the time
2764 taken by this routine is O (n_col), that is, linear in the number of
2765 columns. Not user-callable.
2768 PRIVATE void order_children
2770 /* === Parameters ======================================================= */
2772 Int n_col, /* number of columns of A */
2773 Colamd_Col Col [], /* of size n_col+1 */
2774 Int p [] /* p [0 ... n_col-1] is the column permutation*/
2777 /* === Local variables ================================================== */
2779 Int i ; /* loop counter for all columns */
2780 Int c ; /* column index */
2781 Int parent ; /* index of column's parent */
2782 Int order ; /* column's order */
2784 /* === Order each non-principal column ================================== */
2786 for (i = 0 ; i < n_col ; i++)
2788 /* find an un-ordered non-principal column */
2789 ASSERT (COL_IS_DEAD (i)) ;
2790 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
2793 /* once found, find its principal parent */
2796 parent = Col [parent].shared1.parent ;
2797 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
2799 /* now, order all un-ordered non-principal columns along path */
2800 /* to this parent. collapse tree at the same time */
2802 /* get order of parent */
2803 order = Col [parent].shared2.order ;
2807 ASSERT (Col [c].shared2.order == EMPTY) ;
2809 /* order this column */
2810 Col [c].shared2.order = order++ ;
2812 Col [c].shared1.parent = parent ;
2814 /* get immediate parent of this column */
2815 c = Col [c].shared1.parent ;
2817 /* continue until we hit an ordered column. There are */
2818 /* guarranteed not to be anymore unordered columns */
2819 /* above an ordered column */
2820 } while (Col [c].shared2.order == EMPTY) ;
2822 /* re-order the super_col parent to largest order for this group */
2823 Col [parent].shared2.order = order ;
2827 /* === Generate the permutation ========================================= */
2829 for (c = 0 ; c < n_col ; c++)
2831 p [Col [c].shared2.order] = c ;
2836 /* ========================================================================== */
2837 /* === detect_super_cols ==================================================== */
2838 /* ========================================================================== */
2841 Detects supercolumns by finding matches between columns in the hash buckets.
2842 Check amongst columns in the set A [row_start ... row_start + row_length-1].
2843 The columns under consideration are currently *not* in the degree lists,
2844 and have already been placed in the hash buckets.
2846 The hash bucket for columns whose hash function is equal to h is stored
2849 if head [h] is >= 0, then head [h] contains a degree list, so:
2851 head [h] is the first column in degree bucket h.
2852 Col [head [h]].headhash gives the first column in hash bucket h.
2854 otherwise, the degree list is empty, and:
2856 -(head [h] + 2) is the first column in hash bucket h.
2858 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
2859 column" pointer. Col [c].shared3.hash is used instead as the hash number
2860 for that column. The value of Col [c].shared4.hash_next is the next column
2861 in the same hash bucket.
2863 Assuming no, or "few" hash collisions, the time taken by this routine is
2864 linear in the sum of the sizes (lengths) of each column whose score has
2865 just been computed in the approximate degree computation.
2869 PRIVATE void detect_super_cols
2871 /* === Parameters ======================================================= */
2874 /* these two parameters are only needed when debugging is enabled: */
2875 Int n_col, /* number of columns of A */
2876 Colamd_Row Row [], /* of size n_row+1 */
2879 Colamd_Col Col [], /* of size n_col+1 */
2880 Int A [], /* row indices of A */
2881 Int head [], /* head of degree lists and hash buckets */
2882 Int row_start, /* pointer to set of columns to check */
2883 Int row_length /* number of columns to check */
2886 /* === Local variables ================================================== */
2888 Int hash ; /* hash value for a column */
2889 Int *rp ; /* pointer to a row */
2890 Int c ; /* a column index */
2891 Int super_c ; /* column index of the column to absorb into */
2892 Int *cp1 ; /* column pointer for column super_c */
2893 Int *cp2 ; /* column pointer for column c */
2894 Int length ; /* length of column super_c */
2895 Int prev_c ; /* column preceding c in hash bucket */
2896 Int i ; /* loop counter */
2897 Int *rp_end ; /* pointer to the end of the row */
2898 Int col ; /* a column index in the row to check */
2899 Int head_column ; /* first column in hash bucket or degree list */
2900 Int first_col ; /* first column in hash bucket */
2902 /* === Consider each column in the row ================================== */
2904 rp = &A [row_start] ;
2905 rp_end = rp + row_length ;
2909 if (COL_IS_DEAD (col))
2914 /* get hash number for this column */
2915 hash = Col [col].shared3.hash ;
2916 ASSERT (hash <= n_col) ;
2918 /* === Get the first column in this hash bucket ===================== */
2920 head_column = head [hash] ;
2921 if (head_column > EMPTY)
2923 first_col = Col [head_column].shared3.headhash ;
2927 first_col = - (head_column + 2) ;
2930 /* === Consider each column in the hash bucket ====================== */
2932 for (super_c = first_col ; super_c != EMPTY ;
2933 super_c = Col [super_c].shared4.hash_next)
2935 ASSERT (COL_IS_ALIVE (super_c)) ;
2936 ASSERT (Col [super_c].shared3.hash == hash) ;
2937 length = Col [super_c].length ;
2939 /* prev_c is the column preceding column c in the hash bucket */
2942 /* === Compare super_c with all columns after it ================ */
2944 for (c = Col [super_c].shared4.hash_next ;
2945 c != EMPTY ; c = Col [c].shared4.hash_next)
2947 ASSERT (c != super_c) ;
2948 ASSERT (COL_IS_ALIVE (c)) ;
2949 ASSERT (Col [c].shared3.hash == hash) ;
2951 /* not identical if lengths or scores are different */
2952 if (Col [c].length != length ||
2953 Col [c].shared2.score != Col [super_c].shared2.score)
2959 /* compare the two columns */
2960 cp1 = &A [Col [super_c].start] ;
2961 cp2 = &A [Col [c].start] ;
2963 for (i = 0 ; i < length ; i++)
2965 /* the columns are "clean" (no dead rows) */
2966 ASSERT (ROW_IS_ALIVE (*cp1)) ;
2967 ASSERT (ROW_IS_ALIVE (*cp2)) ;
2968 /* row indices will same order for both supercols, */
2969 /* no gather scatter nessasary */
2970 if (*cp1++ != *cp2++)
2976 /* the two columns are different if the for-loop "broke" */
2983 /* === Got it! two columns are identical =================== */
2985 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
2987 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
2988 Col [c].shared1.parent = super_c ;
2989 KILL_NON_PRINCIPAL_COL (c) ;
2990 /* order c later, in order_children() */
2991 Col [c].shared2.order = EMPTY ;
2992 /* remove c from hash bucket */
2993 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
2997 /* === Empty this hash bucket ======================================= */
2999 if (head_column > EMPTY)
3001 /* corresponding degree list "hash" is not empty */
3002 Col [head_column].shared3.headhash = EMPTY ;
3006 /* corresponding degree list "hash" is empty */
3007 head [hash] = EMPTY ;
3013 /* ========================================================================== */
3014 /* === garbage_collection =================================================== */
3015 /* ========================================================================== */
3018 Defragments and compacts columns and rows in the workspace A. Used when
3019 all avaliable memory has been used while performing row merging. Returns
3020 the index of the first free position in A, after garbage collection. The
3021 time taken by this routine is linear is the size of the array A, which is
3022 itself linear in the number of nonzeros in the input matrix.
3026 PRIVATE Int garbage_collection /* returns the new value of pfree */
3028 /* === Parameters ======================================================= */
3030 Int n_row, /* number of rows */
3031 Int n_col, /* number of columns */
3032 Colamd_Row Row [], /* row info */
3033 Colamd_Col Col [], /* column info */
3034 Int A [], /* A [0 ... Alen-1] holds the matrix */
3035 Int *pfree /* &A [0] ... pfree is in use */
3038 /* === Local variables ================================================== */
3040 Int *psrc ; /* source pointer */
3041 Int *pdest ; /* destination pointer */
3042 Int j ; /* counter */
3043 Int r ; /* a row index */
3044 Int c ; /* a column index */
3045 Int length ; /* length of a row or column */
3049 DEBUG2 (("Defrag..\n")) ;
3050 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
3054 /* === Defragment the columns =========================================== */
3057 for (c = 0 ; c < n_col ; c++)
3059 if (COL_IS_ALIVE (c))
3061 psrc = &A [Col [c].start] ;
3063 /* move and compact the column */
3064 ASSERT (pdest <= psrc) ;
3065 Col [c].start = (Int) (pdest - &A [0]) ;
3066 length = Col [c].length ;
3067 for (j = 0 ; j < length ; j++)
3070 if (ROW_IS_ALIVE (r))
3075 Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
3079 /* === Prepare to defragment the rows =================================== */
3081 for (r = 0 ; r < n_row ; r++)
3083 if (ROW_IS_DEAD (r) || (Row [r].length == 0))
3085 /* This row is already dead, or is of zero length. Cannot compact
3086 * a row of zero length, so kill it. NOTE: in the current version,
3087 * there are no zero-length live rows. Kill the row (for the first
3088 * time, or again) just to be safe. */
3093 /* save first column index in Row [r].shared2.first_column */
3094 psrc = &A [Row [r].start] ;
3095 Row [r].shared2.first_column = *psrc ;
3096 ASSERT (ROW_IS_ALIVE (r)) ;
3097 /* flag the start of the row with the one's complement of row */
3098 *psrc = ONES_COMPLEMENT (r) ;
3105 /* === Defragment the rows ============================================== */
3108 while (psrc < pfree)
3110 /* find a negative number ... the start of a row */
3114 /* get the row index */
3115 r = ONES_COMPLEMENT (*psrc) ;
3116 ASSERT (r >= 0 && r < n_row) ;
3117 /* restore first column index */
3118 *psrc = Row [r].shared2.first_column ;
3119 ASSERT (ROW_IS_ALIVE (r)) ;
3120 ASSERT (Row [r].length > 0) ;
3121 /* move and compact the row */
3122 ASSERT (pdest <= psrc) ;
3123 Row [r].start = (Int) (pdest - &A [0]) ;
3124 length = Row [r].length ;
3125 for (j = 0 ; j < length ; j++)
3128 if (COL_IS_ALIVE (c))
3133 Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
3134 ASSERT (Row [r].length > 0) ;
3140 /* ensure we found all the rows */
3141 ASSERT (debug_rows == 0) ;
3143 /* === Return the new value of pfree ==================================== */
3145 return ((Int) (pdest - &A [0])) ;
3149 /* ========================================================================== */
3150 /* === clear_mark =========================================================== */
3151 /* ========================================================================== */
3154 Clears the Row [].shared2.mark array, and returns the new tag_mark.
3155 Return value is the new tag_mark. Not user-callable.
3158 PRIVATE Int clear_mark /* return the new value for tag_mark */
3160 /* === Parameters ======================================================= */
3162 Int tag_mark, /* new value of tag_mark */
3163 Int max_mark, /* max allowed value of tag_mark */
3165 Int n_row, /* number of rows in A */
3166 Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
3169 /* === Local variables ================================================== */
3173 if (tag_mark <= 0 || tag_mark >= max_mark)
3175 for (r = 0 ; r < n_row ; r++)
3177 if (ROW_IS_ALIVE (r))
3179 Row [r].shared2.mark = 0 ;
3189 /* ========================================================================== */
3190 /* === print_report ========================================================= */
3191 /* ========================================================================== */
3193 PRIVATE void print_report
3196 Int stats [COLAMD_STATS]
3202 PRINTF (("\n%s version %d.%d, %s: ", method,
3203 COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
3207 PRINTF (("No statistics available.\n")) ;
3211 i1 = stats [COLAMD_INFO1] ;
3212 i2 = stats [COLAMD_INFO2] ;
3213 i3 = stats [COLAMD_INFO3] ;
3215 if (stats [COLAMD_STATUS] >= 0)
3221 PRINTF (("ERROR. ")) ;
3224 switch (stats [COLAMD_STATUS])
3227 case COLAMD_OK_BUT_JUMBLED:
3229 PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
3231 PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
3234 PRINTF(("%s: last seen duplicate or out-of-order row index: %d\n",
3235 method, INDEX (i2))) ;
3237 PRINTF(("%s: last seen in column: %d",
3238 method, INDEX (i1))) ;
3240 /* no break - fall through to next case instead */
3246 PRINTF(("%s: number of dense or empty rows ignored: %d\n",
3247 method, stats [COLAMD_DENSE_ROW])) ;
3249 PRINTF(("%s: number of dense or empty columns ignored: %d\n",
3250 method, stats [COLAMD_DENSE_COL])) ;
3252 PRINTF(("%s: number of garbage collections performed: %d\n",
3253 method, stats [COLAMD_DEFRAG_COUNT])) ;
3256 case COLAMD_ERROR_A_not_present:
3258 PRINTF(("Array A (row indices of matrix) not present.\n")) ;
3261 case COLAMD_ERROR_p_not_present:
3263 PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
3266 case COLAMD_ERROR_nrow_negative:
3268 PRINTF(("Invalid number of rows (%d).\n", i1)) ;
3271 case COLAMD_ERROR_ncol_negative:
3273 PRINTF(("Invalid number of columns (%d).\n", i1)) ;
3276 case COLAMD_ERROR_nnz_negative:
3278 PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
3281 case COLAMD_ERROR_p0_nonzero:
3283 PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
3286 case COLAMD_ERROR_A_too_small:
3288 PRINTF(("Array A too small.\n")) ;
3289 PRINTF((" Need Alen >= %d, but given only Alen = %d.\n",
3293 case COLAMD_ERROR_col_length_negative:
3296 (("Column %d has a negative number of nonzero entries (%d).\n",
3300 case COLAMD_ERROR_row_index_out_of_bounds:
3303 (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
3304 INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
3307 case COLAMD_ERROR_out_of_memory:
3309 PRINTF(("Out of memory.\n")) ;
3312 /* v2.4: internal-error case deleted */
3319 /* ========================================================================== */
3320 /* === colamd debugging routines ============================================ */
3321 /* ========================================================================== */
3323 /* When debugging is disabled, the remainder of this file is ignored. */
3328 /* ========================================================================== */
3329 /* === debug_structures ===================================================== */
3330 /* ========================================================================== */
3333 At this point, all empty rows and columns are dead. All live columns
3334 are "clean" (containing no dead rows) and simplicial (no supercolumns
3335 yet). Rows may contain dead columns, but all live rows contain at
3336 least one live column.
3339 PRIVATE void debug_structures
3341 /* === Parameters ======================================================= */
3351 /* === Local variables ================================================== */
3364 /* === Check A, Row, and Col ============================================ */
3366 for (c = 0 ; c < n_col ; c++)
3368 if (COL_IS_ALIVE (c))
3370 len = Col [c].length ;
3371 score = Col [c].shared2.score ;
3372 DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
3374 ASSERT (score >= 0) ;
3375 ASSERT (Col [c].shared1.thickness == 1) ;
3376 cp = &A [Col [c].start] ;
3381 ASSERT (ROW_IS_ALIVE (r)) ;
3386 i = Col [c].shared2.order ;
3387 ASSERT (i >= n_col2 && i < n_col) ;
3391 for (r = 0 ; r < n_row ; r++)
3393 if (ROW_IS_ALIVE (r))
3396 len = Row [r].length ;
3397 deg = Row [r].shared1.degree ;
3400 rp = &A [Row [r].start] ;
3405 if (COL_IS_ALIVE (c))
3416 /* ========================================================================== */
3417 /* === debug_deg_lists ====================================================== */
3418 /* ========================================================================== */
3421 Prints the contents of the degree lists. Counts the number of columns
3422 in the degree list and compares it to the total it should have. Also
3423 checks the row degrees.
3426 PRIVATE void debug_deg_lists
3428 /* === Parameters ======================================================= */
3440 /* === Local variables ================================================== */
3447 /* === Check the degree lists =========================================== */
3449 if (n_col > 10000 && colamd_debug <= 0)
3454 DEBUG4 (("Degree lists: %d\n", min_score)) ;
3455 for (deg = 0 ; deg <= n_col ; deg++)
3462 DEBUG4 (("%d:", deg)) ;
3463 while (col != EMPTY)
3465 DEBUG4 ((" %d", col)) ;
3466 have += Col [col].shared1.thickness ;
3467 ASSERT (COL_IS_ALIVE (col)) ;
3468 col = Col [col].shared4.degree_next ;
3472 DEBUG4 (("should %d have %d\n", should, have)) ;
3473 ASSERT (should == have) ;
3475 /* === Check the row degrees ============================================ */
3477 if (n_row > 10000 && colamd_debug <= 0)
3481 for (row = 0 ; row < n_row ; row++)
3483 if (ROW_IS_ALIVE (row))
3485 ASSERT (Row [row].shared1.degree <= max_deg) ;
3491 /* ========================================================================== */
3492 /* === debug_mark =========================================================== */
3493 /* ========================================================================== */
3496 Ensures that the tag_mark is less that the maximum and also ensures that
3497 each entry in the mark array is less than the tag mark.
3500 PRIVATE void debug_mark
3502 /* === Parameters ======================================================= */
3510 /* === Local variables ================================================== */
3514 /* === Check the Row marks ============================================== */
3516 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
3517 if (n_row > 10000 && colamd_debug <= 0)
3521 for (r = 0 ; r < n_row ; r++)
3523 ASSERT (Row [r].shared2.mark < tag_mark) ;
3528 /* ========================================================================== */
3529 /* === debug_matrix ========================================================= */
3530 /* ========================================================================== */
3533 Prints out the contents of the columns and the rows.
3536 PRIVATE void debug_matrix
3538 /* === Parameters ======================================================= */
3547 /* === Local variables ================================================== */
3556 /* === Dump the rows and columns of the matrix ========================== */
3558 if (colamd_debug < 3)
3562 DEBUG3 (("DUMP MATRIX:\n")) ;
3563 for (r = 0 ; r < n_row ; r++)
3565 DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
3566 if (ROW_IS_DEAD (r))
3570 DEBUG3 (("start %d length %d degree %d\n",
3571 Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
3572 rp = &A [Row [r].start] ;
3573 rp_end = rp + Row [r].length ;
3577 DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ;
3581 for (c = 0 ; c < n_col ; c++)
3583 DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
3584 if (COL_IS_DEAD (c))
3588 DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
3589 Col [c].start, Col [c].length,
3590 Col [c].shared1.thickness, Col [c].shared2.score)) ;
3591 cp = &A [Col [c].start] ;
3592 cp_end = cp + Col [c].length ;
3596 DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ;
3601 PRIVATE void colamd_get_debug
3607 colamd_debug = 0 ; /* no debug printing */
3608 f = fopen ("debug", "r") ;
3609 if (f == (FILE *) NULL)
3615 fscanf (f, "%d", &colamd_debug) ;
3618 DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
3619 method, colamd_debug)) ;