C ALGORITHM 836, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 30, NO. 3, September, 2004, P. 377--380. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # C/ # C/Drivers/ # C/Drivers/colamd_example.c # C/Drivers/colamd_example.out # C/Src/ # C/Src/colamd.c # C/Src/colamd.h # ChangeLog # Makefile # Matlab/ # Matlab/Drivers/ # Matlab/Drivers/colamd_demo.m # Matlab/Drivers/colamd_test.m # Matlab/Drivers/colamdtestmex.c # Matlab/Drivers/luflops.m # Matlab/Drivers/symamdtestmex.c # Matlab/Src/ # Matlab/Src/colamd.m # Matlab/Src/colamdmex.c # Matlab/Src/symamd.m # Matlab/Src/symamdmex.c # README # This archive created: Mon Mar 7 16:05:45 2005 export PATH; PATH=/bin:$PATH if test ! -d 'C' then mkdir 'C' fi cd 'C' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'colamd_example.c' then echo shar: will not over-write existing file "'colamd_example.c'" else cat << "SHAR_EOF" > 'colamd_example.c' /* ========================================================================== */ /* === colamd and symamd example ============================================ */ /* ========================================================================== */ /* colamd example of use, to order the columns of a 5-by-4 matrix with 11 nonzero entries in the following nonzero pattern, with default knobs. x 0 x 0 x 0 x x 0 x x 0 0 0 x x x x 0 0 symamd example of use, to order the rows and columns of a 5-by-5 matrix with 13 nonzero entries in the following nonzero pattern, with default knobs. x x 0 0 0 x x x x 0 0 x x 0 0 0 x 0 x x 0 0 0 x x (where x denotes a nonzero value). September 8, 2003. Version 2.3. See http://www.cise.ufl.edu/research/sparse/colamd/ (the colamd.c file) for the routines this program calls, and for the License. */ /* ========================================================================== */ #include #include "colamd.h" #define A_NNZ 11 #define A_NROW 5 #define A_NCOL 4 #define ALEN (COLAMD_RECOMMENDED (A_NNZ, A_NCOL, A_NROW)) #define B_NNZ 4 #define B_N 5 int main (int argc, char **argv) { /* ====================================================================== */ /* input matrix A definition */ /* ====================================================================== */ int A [ALEN] = { 0, 1, 4, /* row indices of nonzeros in column 0 */ 2, 4, /* row indices of nonzeros in column 1 */ 0, 1, 2, 3, /* row indices of nonzeros in column 2 */ 1, 3} ; /* row indices of nonzeros in column 3 */ int p [ ] = { 0, /* column 0 is in A [0..2] */ 3, /* column 1 is in A [3..4] */ 5, /* column 2 is in A [5..8] */ 9, /* column 3 is in A [9..10] */ A_NNZ} ; /* number of nonzeros in A */ /* ====================================================================== */ /* input matrix B definition */ /* ====================================================================== */ int B [ ] = { /* Note: only strictly lower triangular part */ /* is included, since symamd ignores the */ /* diagonal and upper triangular part of B. */ 1, /* row indices of nonzeros in column 0 */ 2, 3, /* row indices of nonzeros in column 1 */ /* row indices of nonzeros in column 2 (none) */ 4 /* row indices of nonzeros in column 3 */ } ; /* row indices of nonzeros in column 4 (none) */ int q [ ] = { 0, /* column 0 is in B [0] */ 1, /* column 1 is in B [1..2] */ 3, /* column 2 is empty */ 3, /* column 3 is in B [3] */ 4, /* column 4 is empty */ B_NNZ} ; /* number of nonzeros in strictly lower B */ /* ====================================================================== */ /* other variable definitions */ /* ====================================================================== */ int perm [B_N+1] ; /* note the size is N+1 */ int stats [COLAMD_STATS] ; /* for colamd and symamd output statistics */ int row, col, pp, length, ok ; /* ====================================================================== */ /* dump the input matrix A */ /* ====================================================================== */ printf ("colamd %d-by-%d input matrix:\n", A_NROW, A_NCOL) ; for (col = 0 ; col < A_NCOL ; col++) { length = p [col+1] - p [col] ; printf ("Column %d, with %d entries:\n", col, length) ; for (pp = p [col] ; pp < p [col+1] ; pp++) { row = A [pp] ; printf (" row %d\n", row) ; } } /* ====================================================================== */ /* order the matrix. Note that this destroys A and overwrites p */ /* ====================================================================== */ ok = colamd (A_NROW, A_NCOL, ALEN, A, p, (double *) NULL, stats) ; colamd_report (stats) ; if (!ok) { printf ("colamd error!\n") ; exit (1) ; } /* ====================================================================== */ /* print the column ordering */ /* ====================================================================== */ printf ("colamd column ordering:\n") ; printf ("1st column: %d\n", p [0]) ; printf ("2nd column: %d\n", p [1]) ; printf ("3rd column: %d\n", p [2]) ; printf ("4th column: %d\n", p [3]) ; /* ====================================================================== */ /* dump the strictly lower triangular part of symmetric input matrix B */ /* ====================================================================== */ printf ("\n\nsymamd %d-by-%d input matrix:\n", B_N, B_N) ; printf ("Entries in strictly lower triangular part:\n") ; for (col = 0 ; col < B_N ; col++) { length = q [col+1] - q [col] ; printf ("Column %d, with %d entries:\n", col, length) ; for (pp = q [col] ; pp < q [col+1] ; pp++) { row = B [pp] ; printf (" row %d\n", row) ; } } /* ====================================================================== */ /* order the matrix B. Note that this does not modify B or q. */ /* ====================================================================== */ ok = symamd (B_N, B, q, perm, (double *) NULL, stats, &calloc, &free) ; symamd_report (stats) ; if (!ok) { printf ("symamd error!\n") ; exit (1) ; } /* ====================================================================== */ /* print the symmetric ordering */ /* ====================================================================== */ printf ("symamd column ordering:\n") ; printf ("1st row/column: %d\n", perm [0]) ; printf ("2nd row/column: %d\n", perm [1]) ; printf ("3rd row/column: %d\n", perm [2]) ; printf ("4th row/column: %d\n", perm [3]) ; printf ("5th row/column: %d\n", perm [4]) ; exit (0) ; } SHAR_EOF fi # end of overwriting check if test -f 'colamd_example.out' then echo shar: will not over-write existing file "'colamd_example.out'" else cat << "SHAR_EOF" > 'colamd_example.out' colamd 5-by-4 input matrix: Column 0, with 3 entries: row 0 row 1 row 4 Column 1, with 2 entries: row 2 row 4 Column 2, with 4 entries: row 0 row 1 row 2 row 3 Column 3, with 2 entries: row 1 row 3 colamd: OK. colamd: number of dense or empty rows ignored: 1 colamd: number of dense or empty columns ignored: 2 colamd: number of garbage collections performed: 0 colamd column ordering: 1st column: 1 2nd column: 3 3rd column: 0 4th column: 2 symamd 5-by-5 input matrix: Entries in strictly lower triangular part: Column 0, with 1 entries: row 1 Column 1, with 2 entries: row 2 row 3 Column 2, with 0 entries: Column 3, with 1 entries: row 4 Column 4, with 0 entries: symamd: OK. symamd: number of dense or empty rows ignored: 0 symamd: number of dense or empty columns ignored: 0 symamd: number of garbage collections performed: 0 symamd column ordering: 1st row/column: 0 2nd row/column: 2 3rd row/column: 1 4th row/column: 3 5th row/column: 4 SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'colamd.c' then echo shar: will not over-write existing file "'colamd.c'" else cat << "SHAR_EOF" > 'colamd.c' /* ========================================================================== */ /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */ /* ========================================================================== */ /* colamd: an approximate minimum degree column ordering algorithm, for LU factorization of symmetric or unsymmetric matrices, QR factorization, least squares, interior point methods for linear programming problems, and other related problems. symamd: an approximate minimum degree ordering algorithm for Cholesky factorization of symmetric matrices. Purpose: Colamd computes a permutation Q such that the Cholesky factorization of (AQ)'(AQ) has less fill-in and requires fewer floating point operations than A'A. This also provides a good ordering for sparse partial pivoting methods, P(AQ) = LU, where Q is computed prior to numerical factorization, and P is computed during numerical factorization via conventional partial pivoting with row interchanges. Colamd is the column ordering method used in SuperLU, part of the ScaLAPACK library. It is also available as built-in function in MATLAB Version 6, available from MathWorks, Inc. (http://www.mathworks.com). This routine can be used in place of colmmd in MATLAB. Symamd computes a permutation P of a symmetric matrix A such that the Cholesky factorization of PAP' has less fill-in and requires fewer floating point operations than A. Symamd constructs a matrix M such that M'M has the same nonzero pattern of A, and then orders the columns of M using colmmd. The column ordering of M is then returned as the row and column ordering P of A. Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Copyright and License: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. Permission is hereby granted to use, copy, modify, and/or distribute this program, provided that the Copyright, this License, and the Availability of the original version is retained on all copies and made accessible to the end-user of any code or package that includes COLAMD or any modified version of COLAMD. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c file. It requires the colamd.h file. It is required by the colamdmex.c and symamdmex.c files, for the MATLAB interface to colamd and symamd. See the ChangeLog file for changes since Version 1.0. */ /* ========================================================================== */ /* === Description of user-callable routines ================================ */ /* ========================================================================== */ /* ---------------------------------------------------------------------------- colamd_recommended: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" int colamd_recommended (int nnz, int n_row, int n_col) ; or as a C macro #include "colamd.h" Alen = COLAMD_RECOMMENDED (int nnz, int n_row, int n_col) ; Purpose: Returns recommended value of Alen for use by colamd. Returns -1 if any input argument is negative. The use of this routine or macro is optional. Note that the macro uses its arguments more than once, so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED. Not needed for symamd, which dynamically allocates its own memory. Arguments (all input arguments): int nnz ; Number of nonzeros in the matrix A. This must be the same value as p [n_col] in the call to colamd - otherwise you will get a wrong value of the recommended memory to use. int n_row ; Number of rows in the matrix A. int n_col ; Number of columns in the matrix A. ---------------------------------------------------------------------------- colamd_set_defaults: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" colamd_set_defaults (double knobs [COLAMD_KNOBS]) ; Purpose: Sets the default parameters. The use of this routine is optional. Arguments: double knobs [COLAMD_KNOBS] ; Output only. Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) entries are removed prior to ordering. Columns with more than (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to ordering, and placed last in the output column ordering. Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0]. Rows and columns with more than (knobs [COLAMD_DENSE_ROW] * n) entries are removed prior to ordering, and placed last in the output ordering. COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1, respectively, in colamd.h. Default values of these two knobs are both 0.5. Currently, only knobs [0] and knobs [1] are used, but future versions may use more knobs. If so, they will be properly set to their defaults by the future version of colamd_set_defaults, so that the code that calls colamd will not need to change, assuming that you either use colamd_set_defaults, or pass a (double *) NULL pointer as the knobs array to colamd or symamd. ---------------------------------------------------------------------------- colamd: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" int colamd (int n_row, int n_col, int Alen, int *A, int *p, double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ; Purpose: Computes a column ordering (Q) of A such that P(AQ)=LU or (AQ)'AQ=LL' have less fill-in and require fewer floating point operations than factorizing the unpermuted matrix A or A'A, respectively. Returns: TRUE (1) if successful, FALSE (0) otherwise. Arguments: int n_row ; Input argument. Number of rows in the matrix A. Restriction: n_row >= 0. Colamd returns FALSE if n_row is negative. int n_col ; Input argument. Number of columns in the matrix A. Restriction: n_col >= 0. Colamd returns FALSE if n_col is negative. int Alen ; Input argument. Restriction (see note): Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col Colamd returns FALSE if these conditions are not met. Note: this restriction makes an modest assumption regarding the size of the two typedef's structures in colamd.h. We do, however, guarantee that Alen >= colamd_recommended (nnz, n_row, n_col) or equivalently as a C preprocessor macro: Alen >= COLAMD_RECOMMENDED (nnz, n_row, n_col) will be sufficient. int A [Alen] ; Input argument, undefined on output. A is an integer array of size Alen. Alen must be at least as large as the bare minimum value given above, but this is very low, and can result in excessive run time. For best performance, we recommend that Alen be greater than or equal to colamd_recommended (nnz, n_row, n_col), which adds nnz/5 to the bare minimum value given above. On input, the row indices of the entries in column c of the matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a given column c need not be in ascending order, and duplicate row indices may be be present. However, colamd will work a little faster if both of these conditions are met (Colamd puts the matrix into this format, if it finds that the the conditions are not met). The matrix is 0-based. That is, rows are in the range 0 to n_row-1, and columns are in the range 0 to n_col-1. Colamd returns FALSE if any row index is out of range. The contents of A are modified during ordering, and are undefined on output. int p [n_col+1] ; Both input and output argument. p is an integer array of size n_col+1. On input, it holds the "pointers" for the column form of the matrix A. Column c of the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first entry, p [0], must be zero, and p [c] <= p [c+1] must hold for all c in the range 0 to n_col-1. The value p [n_col] is thus the total number of entries in the pattern of the matrix A. Colamd returns FALSE if these conditions are not met. On output, if colamd returns TRUE, the array p holds the column permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is the first column index in the new ordering, and p [n_col-1] is the last. That is, p [k] = j means that column j of A is the kth pivot column, in AQ, where k is in the range 0 to n_col-1 (p [0] = j means that column j of A is the first column in AQ). If colamd returns FALSE, then no permutation is returned, and p is undefined on output. double knobs [COLAMD_KNOBS] ; Input argument. See colamd_set_defaults for a description. int stats [COLAMD_STATS] ; Output argument. Statistics on the ordering, and error status. See colamd.h for related definitions. Colamd returns FALSE if stats is not present. stats [0]: number of dense or empty rows ignored. stats [1]: number of dense or empty columns ignored (and ordered last in the output permutation p) Note that a row can become "empty" if it contains only "dense" and/or "empty" columns, and similarly a column can become "empty" if it only contains "dense" and/or "empty" rows. stats [2]: number of garbage collections performed. This can be excessively high if Alen is close to the minimum required value. stats [3]: status code. < 0 is an error code. > 1 is a warning or notice. 0 OK. Each column of the input matrix contained row indices in increasing order, with no duplicates. 1 OK, but columns of input matrix were jumbled (unsorted columns or duplicate entries). Colamd had to do some extra work to sort the matrix first and remove duplicate entries, but it still was able to return a valid permutation (return value of colamd was TRUE). stats [4]: highest numbered column that is unsorted or has duplicate entries. stats [5]: last seen duplicate or unsorted row index. stats [6]: number of duplicate or unsorted row indices. -1 A is a null pointer -2 p is a null pointer -3 n_row is negative stats [4]: n_row -4 n_col is negative stats [4]: n_col -5 number of nonzeros in matrix is negative stats [4]: number of nonzeros, p [n_col] -6 p [0] is nonzero stats [4]: p [0] -7 A is too small stats [4]: required size stats [5]: actual size (Alen) -8 a column has a negative number of entries stats [4]: column with < 0 entries stats [5]: number of entries in col -9 a row index is out of bounds stats [4]: column with bad row index stats [5]: bad row index stats [6]: n_row, # of rows of matrx -10 (unused; see symamd.c) -999 (unused; see symamd.c) Future versions may return more statistics in the stats array. Example: See http://www.cise.ufl.edu/research/sparse/colamd/example.c for a complete example. To order the columns of a 5-by-4 matrix with 11 nonzero entries in the following nonzero pattern x 0 x 0 x 0 x x 0 x x 0 0 0 x x x x 0 0 with default knobs and no output statistics, do the following: #include "colamd.h" #define ALEN COLAMD_RECOMMENDED (11, 5, 4) int A [ALEN] = {1, 2, 5, 3, 5, 1, 2, 3, 4, 2, 4} ; int p [ ] = {0, 3, 5, 9, 11} ; int stats [COLAMD_STATS] ; colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ; The permutation is returned in the array p, and A is destroyed. ---------------------------------------------------------------------------- symamd: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" int symamd (int n, int *A, int *p, int *perm, double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS], void (*allocate) (size_t, size_t), void (*release) (void *)) ; Purpose: The symamd routine computes an ordering P of a symmetric sparse matrix A such that the Cholesky factorization PAP' = LL' remains sparse. It is based on a column ordering of a matrix M constructed so that the nonzero pattern of M'M is the same as A. The matrix A is assumed to be symmetric; only the strictly lower triangular part is accessed. You must pass your selected memory allocator (usually calloc/free or mxCalloc/mxFree) to symamd, for it to allocate memory for the temporary matrix M. Returns: TRUE (1) if successful, FALSE (0) otherwise. Arguments: int n ; Input argument. Number of rows and columns in the symmetrix matrix A. Restriction: n >= 0. Symamd returns FALSE if n is negative. int A [nnz] ; Input argument. A is an integer array of size nnz, where nnz = p [n]. The row indices of the entries in column c of the matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a given column c need not be in ascending order, and duplicate row indices may be present. However, symamd will run faster if the columns are in sorted order with no duplicate entries. The matrix is 0-based. That is, rows are in the range 0 to n-1, and columns are in the range 0 to n-1. Symamd returns FALSE if any row index is out of range. The contents of A are not modified. int p [n+1] ; Input argument. p is an integer array of size n+1. On input, it holds the "pointers" for the column form of the matrix A. Column c of the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first entry, p [0], must be zero, and p [c] <= p [c+1] must hold for all c in the range 0 to n-1. The value p [n] is thus the total number of entries in the pattern of the matrix A. Symamd returns FALSE if these conditions are not met. The contents of p are not modified. int perm [n+1] ; Output argument. On output, if symamd returns TRUE, the array perm holds the permutation P, where perm [0] is the first index in the new ordering, and perm [n-1] is the last. That is, perm [k] = j means that row and column j of A is the kth column in PAP', where k is in the range 0 to n-1 (perm [0] = j means that row and column j of A are the first row and column in PAP'). The array is used as a workspace during the ordering, which is why it must be of length n+1, not just n. double knobs [COLAMD_KNOBS] ; Input argument. See colamd_set_defaults for a description. int stats [COLAMD_STATS] ; Output argument. Statistics on the ordering, and error status. See colamd.h for related definitions. Symamd returns FALSE if stats is not present. stats [0]: number of dense or empty row and columns ignored (and ordered last in the output permutation perm). Note that a row/column can become "empty" if it contains only "dense" and/or "empty" columns/rows. stats [1]: (same as stats [0]) stats [2]: number of garbage collections performed. stats [3]: status code. < 0 is an error code. > 1 is a warning or notice. 0 OK. Each column of the input matrix contained row indices in increasing order, with no duplicates. 1 OK, but columns of input matrix were jumbled (unsorted columns or duplicate entries). Symamd had to do some extra work to sort the matrix first and remove duplicate entries, but it still was able to return a valid permutation (return value of symamd was TRUE). stats [4]: highest numbered column that is unsorted or has duplicate entries. stats [5]: last seen duplicate or unsorted row index. stats [6]: number of duplicate or unsorted row indices. -1 A is a null pointer -2 p is a null pointer -3 (unused, see colamd.c) -4 n is negative stats [4]: n -5 number of nonzeros in matrix is negative stats [4]: # of nonzeros (p [n]). -6 p [0] is nonzero stats [4]: p [0] -7 (unused) -8 a column has a negative number of entries stats [4]: column with < 0 entries stats [5]: number of entries in col -9 a row index is out of bounds stats [4]: column with bad row index stats [5]: bad row index stats [6]: n_row, # of rows of matrx -10 out of memory (unable to allocate temporary workspace for M or count arrays using the "allocate" routine passed into symamd). -999 internal error. colamd failed to order the matrix M, when it should have succeeded. This indicates a bug. If this (and *only* this) error code occurs, please contact the authors. Don't contact the authors if you get any other error code. Future versions may return more statistics in the stats array. void * (*allocate) (size_t, size_t) A pointer to a function providing memory allocation. The allocated memory must be returned initialized to zero. For a C application, this argument should normally be a pointer to calloc. For a MATLAB mexFunction, the routine mxCalloc is passed instead. void (*release) (size_t, size_t) A pointer to a function that frees memory allocated by the memory allocation routine above. For a C application, this argument should normally be a pointer to free. For a MATLAB mexFunction, the routine mxFree is passed instead. ---------------------------------------------------------------------------- colamd_report: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" colamd_report (int stats [COLAMD_STATS]) ; Purpose: Prints the error status and statistics recorded in the stats array on the standard error output (for a standard C routine) or on the MATLAB output (for a mexFunction). Arguments: int stats [COLAMD_STATS] ; Input only. Statistics from colamd. ---------------------------------------------------------------------------- symamd_report: ---------------------------------------------------------------------------- C syntax: #include "colamd.h" symamd_report (int stats [COLAMD_STATS]) ; Purpose: Prints the error status and statistics recorded in the stats array on the standard error output (for a standard C routine) or on the MATLAB output (for a mexFunction). Arguments: int stats [COLAMD_STATS] ; Input only. Statistics from symamd. */ /* ========================================================================== */ /* === Scaffolding code definitions ======================================== */ /* ========================================================================== */ /* Ensure that debugging is turned off: */ #ifndef NDEBUG #define NDEBUG #endif /* NDEBUG */ /* Our "scaffolding code" philosophy: In our opinion, well-written library code should keep its "debugging" code, and just normally have it turned off by the compiler so as not to interfere with performance. This serves several purposes: (1) assertions act as comments to the reader, telling you what the code expects at that point. All assertions will always be true (unless there really is a bug, of course). (2) leaving in the scaffolding code assists anyone who would like to modify the code, or understand the algorithm (by reading the debugging output, one can get a glimpse into what the code is doing). (3) (gasp!) for actually finding bugs. This code has been heavily tested and "should" be fully functional and bug-free ... but you never know... To enable debugging, comment out the "#define NDEBUG" above. For a MATLAB mexFunction, you will also need to modify mexopts.sh to remove the -DNDEBUG definition. The code will become outrageously slow when debugging is enabled. To control the level of debugging output, set an environment variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging, you should see the following message on the standard output: colamd: debug version, D = 1 (THIS WILL BE SLOW!) or a similar message for symamd. If you don't, then debugging has not been enabled. */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "colamd.h" #include #ifdef MATLAB_MEX_FILE #include "mex.h" #include "matrix.h" #else #include #include #endif /* MATLAB_MEX_FILE */ /* ========================================================================== */ /* === Definitions ========================================================== */ /* ========================================================================== */ /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */ #define PUBLIC #define PRIVATE static #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #define MIN(a,b) (((a) < (b)) ? (a) : (b)) #define ONES_COMPLEMENT(r) (-(r)-1) /* -------------------------------------------------------------------------- */ /* Change for version 2.1: define TRUE and FALSE only if not yet defined */ /* -------------------------------------------------------------------------- */ #ifndef TRUE #define TRUE (1) #endif #ifndef FALSE #define FALSE (0) #endif /* -------------------------------------------------------------------------- */ #define EMPTY (-1) /* Row and column status */ #define ALIVE (0) #define DEAD (-1) /* Column status */ #define DEAD_PRINCIPAL (-1) #define DEAD_NON_PRINCIPAL (-2) /* Macros for row and column status update and checking. */ #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) #define COL_IS_DEAD(c) (Col [c].start < ALIVE) #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } /* ========================================================================== */ /* === Colamd reporting mechanism =========================================== */ /* ========================================================================== */ #ifdef MATLAB_MEX_FILE /* use mexPrintf in a MATLAB mexFunction, for debugging and statistics output */ #define PRINTF mexPrintf /* In MATLAB, matrices are 1-based to the user, but 0-based internally */ #define INDEX(i) ((i)+1) #else /* Use printf in standard C environment, for debugging and statistics output. */ /* Output is generated only if debugging is enabled at compile time, or if */ /* the caller explicitly calls colamd_report or symamd_report. */ #define PRINTF printf /* In C, matrices are 0-based and indices are reported as such in *_report */ #define INDEX(i) (i) #endif /* MATLAB_MEX_FILE */ /* ========================================================================== */ /* === Prototypes of PRIVATE routines ======================================= */ /* ========================================================================== */ PRIVATE int init_rows_cols ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [], int p [], int stats [COLAMD_STATS] ) ; PRIVATE void init_scoring ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [], int head [], double knobs [COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg ) ; PRIVATE int find_ordering ( int n_row, int n_col, int Alen, Colamd_Row Row [], Colamd_Col Col [], int A [], int head [], int n_col2, int max_deg, int pfree ) ; PRIVATE void order_children ( int n_col, Colamd_Col Col [], int p [] ) ; PRIVATE void detect_super_cols ( #ifndef NDEBUG int n_col, Colamd_Row Row [], #endif /* NDEBUG */ Colamd_Col Col [], int A [], int head [], int row_start, int row_length ) ; PRIVATE int garbage_collection ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [], int *pfree ) ; PRIVATE int clear_mark ( int n_row, Colamd_Row Row [] ) ; PRIVATE void print_report ( char *method, int stats [COLAMD_STATS] ) ; /* ========================================================================== */ /* === Debugging prototypes and definitions ================================= */ /* ========================================================================== */ #ifndef NDEBUG /* colamd_debug is the *ONLY* global variable, and is only */ /* present when debugging */ PRIVATE int colamd_debug ; /* debug print level */ #define DEBUG0(params) { (void) PRINTF params ; } #define DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; } #define DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; } #define DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; } #define DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; } #ifdef MATLAB_MEX_FILE #define ASSERT(expression) (mxAssert ((expression), "")) #else #define ASSERT(expression) (assert (expression)) #endif /* MATLAB_MEX_FILE */ PRIVATE void colamd_get_debug /* gets the debug print level from getenv */ ( char *method ) ; PRIVATE void debug_deg_lists ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int head [], int min_score, int should, int max_deg ) ; PRIVATE void debug_mark ( int n_row, Colamd_Row Row [], int tag_mark, int max_mark ) ; PRIVATE void debug_matrix ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [] ) ; PRIVATE void debug_structures ( int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [], int n_col2 ) ; #else /* NDEBUG */ /* === No debugging ========================================================= */ #define DEBUG0(params) ; #define DEBUG1(params) ; #define DEBUG2(params) ; #define DEBUG3(params) ; #define DEBUG4(params) ; #define ASSERT(expression) ((void) 0) #endif /* NDEBUG */ /* ========================================================================== */ /* ========================================================================== */ /* === USER-CALLABLE ROUTINES: ============================================== */ /* ========================================================================== */ /* ========================================================================== */ /* === colamd_recommended =================================================== */ /* ========================================================================== */ /* The colamd_recommended routine returns the suggested size for Alen. This value has been determined to provide good balance between the number of garbage collections and the memory requirements for colamd. If any argument is negative, a -1 is returned as an error condition. This function is also available as a macro defined in colamd.h, so that you can use it for a statically-allocated array size. */ PUBLIC int colamd_recommended /* returns recommended value of Alen. */ ( /* === Parameters ======================================================= */ int nnz, /* number of nonzeros in A */ int n_row, /* number of rows in A */ int n_col /* number of columns in A */ ) { return (COLAMD_RECOMMENDED (nnz, n_row, n_col)) ; } /* ========================================================================== */ /* === colamd_set_defaults ================================================== */ /* ========================================================================== */ /* The colamd_set_defaults routine sets the default values of the user- controllable parameters for colamd: knobs [0] rows with knobs[0]*n_col entries or more are removed prior to ordering in colamd. Rows and columns with knobs[0]*n_col entries or more are removed prior to ordering in symamd and placed last in the output ordering. knobs [1] columns with knobs[1]*n_row entries or more are removed prior to ordering in colamd, and placed last in the column permutation. Symamd ignores this knob. knobs [2..19] unused, but future versions might use this */ PUBLIC void colamd_set_defaults ( /* === Parameters ======================================================= */ double knobs [COLAMD_KNOBS] /* knob array */ ) { /* === Local variables ================================================== */ int i ; if (!knobs) { return ; /* no knobs to initialize */ } for (i = 0 ; i < COLAMD_KNOBS ; i++) { knobs [i] = 0 ; } knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ } /* ========================================================================== */ /* === symamd =============================================================== */ /* ========================================================================== */ PUBLIC int symamd /* return TRUE if OK, FALSE otherwise */ ( /* === Parameters ======================================================= */ int n, /* number of rows and columns of A */ int A [], /* row indices of A */ int p [], /* column pointers of A */ int perm [], /* output permutation, size n+1 */ double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */ int stats [COLAMD_STATS], /* output statistics and error codes */ void * (*allocate) (size_t, size_t), /* pointer to calloc (ANSI C) or */ /* mxCalloc (for MATLAB mexFunction) */ void (*release) (void *) /* pointer to free (ANSI C) or */ /* mxFree (for MATLAB mexFunction) */ ) { /* === Local variables ================================================== */ int *count ; /* length of each column of M, and col pointer*/ int *mark ; /* mark array for finding duplicate entries */ int *M ; /* row indices of matrix M */ int Mlen ; /* length of M */ int n_row ; /* number of rows in M */ int nnz ; /* number of entries in A */ int i ; /* row index of A */ int j ; /* column index of A */ int k ; /* row index of M */ int mnz ; /* number of nonzeros in M */ int pp ; /* index into a column of A */ int last_row ; /* last row seen in the current column */ int length ; /* number of nonzeros in a column */ double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */ double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */ int cstats [COLAMD_STATS] ; /* colamd stats */ #ifndef NDEBUG colamd_get_debug ("symamd") ; #endif /* NDEBUG */ /* === Check the input arguments ======================================== */ if (!stats) { DEBUG0 (("symamd: stats not present\n")) ; return (FALSE) ; } for (i = 0 ; i < COLAMD_STATS ; i++) { stats [i] = 0 ; } stats [COLAMD_STATUS] = COLAMD_OK ; stats [COLAMD_INFO1] = -1 ; stats [COLAMD_INFO2] = -1 ; if (!A) { stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; DEBUG0 (("symamd: A not present\n")) ; return (FALSE) ; } if (!p) /* p is not present */ { stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; DEBUG0 (("symamd: p not present\n")) ; return (FALSE) ; } if (n < 0) /* n must be >= 0 */ { stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; stats [COLAMD_INFO1] = n ; DEBUG0 (("symamd: n negative %d\n", n)) ; return (FALSE) ; } nnz = p [n] ; if (nnz < 0) /* nnz must be >= 0 */ { stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; stats [COLAMD_INFO1] = nnz ; DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ; return (FALSE) ; } if (p [0] != 0) { stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; stats [COLAMD_INFO1] = p [0] ; DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ; return (FALSE) ; } /* === If no knobs, set default knobs =================================== */ if (!knobs) { colamd_set_defaults (default_knobs) ; knobs = default_knobs ; } /* === Allocate count and mark ========================================== */ count = (int *) ((*allocate) (n+1, sizeof (int))) ; if (!count) { stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ; return (FALSE) ; } mark = (int *) ((*allocate) (n+1, sizeof (int))) ; if (!mark) { stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; (*release) ((void *) count) ; DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ; return (FALSE) ; } /* === Compute column counts of M, check if A is valid ================== */ stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ for (i = 0 ; i < n ; i++) { mark [i] = -1 ; } for (j = 0 ; j < n ; j++) { last_row = -1 ; length = p [j+1] - p [j] ; if (length < 0) { /* column pointers must be non-decreasing */ stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; stats [COLAMD_INFO1] = j ; stats [COLAMD_INFO2] = length ; (*release) ((void *) count) ; (*release) ((void *) mark) ; DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ; return (FALSE) ; } for (pp = p [j] ; pp < p [j+1] ; pp++) { i = A [pp] ; if (i < 0 || i >= n) { /* row index i, in column j, is out of bounds */ stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; stats [COLAMD_INFO1] = j ; stats [COLAMD_INFO2] = i ; stats [COLAMD_INFO3] = n ; (*release) ((void *) count) ; (*release) ((void *) mark) ; DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ; return (FALSE) ; } if (i <= last_row || mark [i] == j) { /* row index is unsorted or repeated (or both), thus col */ /* is jumbled. This is a notice, not an error condition. */ stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; stats [COLAMD_INFO1] = j ; stats [COLAMD_INFO2] = i ; (stats [COLAMD_INFO3]) ++ ; DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ; } if (i > j && mark [i] != j) { /* row k of M will contain column indices i and j */ count [i]++ ; count [j]++ ; } /* mark the row as having been seen in this column */ mark [i] = j ; last_row = i ; } } if (stats [COLAMD_STATUS] == COLAMD_OK) { /* if there are no duplicate entries, then mark is no longer needed */ (*release) ((void *) mark) ; } /* === Compute column pointers of M ===================================== */ /* use output permutation, perm, for column pointers of M */ perm [0] = 0 ; for (j = 1 ; j <= n ; j++) { perm [j] = perm [j-1] + count [j-1] ; } for (j = 0 ; j < n ; j++) { count [j] = perm [j] ; } /* === Construct M ====================================================== */ mnz = perm [n] ; n_row = mnz / 2 ; Mlen = colamd_recommended (mnz, n_row, n) ; M = (int *) ((*allocate) (Mlen, sizeof (int))) ; DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %d\n", n_row, n, mnz, Mlen)) ; if (!M) { stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; (*release) ((void *) count) ; (*release) ((void *) mark) ; DEBUG0 (("symamd: allocate M (size %d) failed\n", Mlen)) ; return (FALSE) ; } k = 0 ; if (stats [COLAMD_STATUS] == COLAMD_OK) { /* Matrix is OK */ for (j = 0 ; j < n ; j++) { ASSERT (p [j+1] - p [j] >= 0) ; for (pp = p [j] ; pp < p [j+1] ; pp++) { i = A [pp] ; ASSERT (i >= 0 && i < n) ; if (i > j) { /* row k of M contains column indices i and j */ M [count [i]++] = k ; M [count [j]++] = k ; k++ ; } } } } else { /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */ DEBUG0 (("symamd: Duplicates in A.\n")) ; for (i = 0 ; i < n ; i++) { mark [i] = -1 ; } for (j = 0 ; j < n ; j++) { ASSERT (p [j+1] - p [j] >= 0) ; for (pp = p [j] ; pp < p [j+1] ; pp++) { i = A [pp] ; ASSERT (i >= 0 && i < n) ; if (i > j && mark [i] != j) { /* row k of M contains column indices i and j */ M [count [i]++] = k ; M [count [j]++] = k ; k++ ; mark [i] = j ; } } } (*release) ((void *) mark) ; } /* count and mark no longer needed */ (*release) ((void *) count) ; ASSERT (k == n_row) ; /* === Adjust the knobs for M =========================================== */ for (i = 0 ; i < COLAMD_KNOBS ; i++) { cknobs [i] = knobs [i] ; } /* there are no dense rows in M */ cknobs [COLAMD_DENSE_ROW] = 1.0 ; if (n_row != 0 && n < n_row) { /* On input, the knob is a fraction of 1..n, the number of rows of A. */ /* Convert it to a fraction of 1..n_row, of the number of rows of M. */ cknobs [COLAMD_DENSE_COL] = (knobs [COLAMD_DENSE_ROW] * n) / n_row ; } else { /* no dense columns in M */ cknobs [COLAMD_DENSE_COL] = 1.0 ; } DEBUG0 (("symamd: dense col knob for M: %g\n", cknobs [COLAMD_DENSE_COL])) ; /* === Order the columns of M =========================================== */ if (!colamd (n_row, n, Mlen, M, perm, cknobs, cstats)) { /* This "cannot" happen, unless there is a bug in the code. */ stats [COLAMD_STATUS] = COLAMD_ERROR_internal_error ; (*release) ((void *) M) ; DEBUG0 (("symamd: internal error!\n")) ; return (FALSE) ; } /* Note that the output permutation is now in perm */ /* === get the statistics for symamd from colamd ======================== */ /* note that a dense column in colamd means a dense row and col in symamd */ stats [COLAMD_DENSE_ROW] = cstats [COLAMD_DENSE_COL] ; stats [COLAMD_DENSE_COL] = cstats [COLAMD_DENSE_COL] ; stats [COLAMD_DEFRAG_COUNT] = cstats [COLAMD_DEFRAG_COUNT] ; /* === Free M =========================================================== */ (*release) ((void *) M) ; DEBUG0 (("symamd: done.\n")) ; return (TRUE) ; } /* ========================================================================== */ /* === colamd =============================================================== */ /* ========================================================================== */ /* The colamd routine computes a column ordering Q of a sparse matrix A such that the LU factorization P(AQ) = LU remains sparse, where P is selected via partial pivoting. The routine can also be viewed as providing a permutation Q such that the Cholesky factorization (AQ)'(AQ) = LL' remains sparse. */ PUBLIC int colamd /* returns TRUE if successful, FALSE otherwise*/ ( /* === Parameters ======================================================= */ int n_row, /* number of rows in A */ int n_col, /* number of columns in A */ int Alen, /* length of A */ int A [], /* row indices of A */ int p [], /* pointers to columns in A */ double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */ int stats [COLAMD_STATS] /* output statistics and error codes */ ) { /* === Local variables ================================================== */ int i ; /* loop index */ int nnz ; /* nonzeros in A */ int Row_size ; /* size of Row [], in integers */ int Col_size ; /* size of Col [], in integers */ int need ; /* minimum required length of A */ Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ int n_col2 ; /* number of non-dense, non-empty columns */ int n_row2 ; /* number of non-dense, non-empty rows */ int ngarbage ; /* number of garbage collections performed */ int max_deg ; /* maximum row degree */ double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ #ifndef NDEBUG colamd_get_debug ("colamd") ; #endif /* NDEBUG */ /* === Check the input arguments ======================================== */ if (!stats) { DEBUG0 (("colamd: stats not present\n")) ; return (FALSE) ; } for (i = 0 ; i < COLAMD_STATS ; i++) { stats [i] = 0 ; } stats [COLAMD_STATUS] = COLAMD_OK ; stats [COLAMD_INFO1] = -1 ; stats [COLAMD_INFO2] = -1 ; if (!A) /* A is not present */ { stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; DEBUG0 (("colamd: A not present\n")) ; return (FALSE) ; } if (!p) /* p is not present */ { stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; DEBUG0 (("colamd: p not present\n")) ; return (FALSE) ; } if (n_row < 0) /* n_row must be >= 0 */ { stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; stats [COLAMD_INFO1] = n_row ; DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; return (FALSE) ; } if (n_col < 0) /* n_col must be >= 0 */ { stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; stats [COLAMD_INFO1] = n_col ; DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; return (FALSE) ; } nnz = p [n_col] ; if (nnz < 0) /* nnz must be >= 0 */ { stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; stats [COLAMD_INFO1] = nnz ; DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; return (FALSE) ; } if (p [0] != 0) { stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; stats [COLAMD_INFO1] = p [0] ; DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; return (FALSE) ; } /* === If no knobs, set default knobs =================================== */ if (!knobs) { colamd_set_defaults (default_knobs) ; knobs = default_knobs ; } /* === Allocate the Row and Col arrays from array A ===================== */ Col_size = COLAMD_C (n_col) ; Row_size = COLAMD_R (n_row) ; need = 2*nnz + n_col + Col_size + Row_size ; if (need > Alen) { /* not enough space in array A to perform the ordering */ stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; stats [COLAMD_INFO1] = need ; stats [COLAMD_INFO2] = Alen ; DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); return (FALSE) ; } Alen -= Col_size + Row_size ; Col = (Colamd_Col *) &A [Alen] ; Row = (Colamd_Row *) &A [Alen + Col_size] ; /* === Construct the row and column data structures ===================== */ if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) { /* input matrix is invalid */ DEBUG0 (("colamd: Matrix invalid\n")) ; return (FALSE) ; } /* === Initialize scores, kill dense rows/columns ======================= */ init_scoring (n_row, n_col, Row, Col, A, p, knobs, &n_row2, &n_col2, &max_deg) ; /* === Order the supercolumns =========================================== */ ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, n_col2, max_deg, 2*nnz) ; /* === Order the non-principal columns ================================== */ order_children (n_col, Col, p) ; /* === Return statistics in stats ======================================= */ stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; stats [COLAMD_DENSE_COL] = n_col - n_col2 ; stats [COLAMD_DEFRAG_COUNT] = ngarbage ; DEBUG0 (("colamd: done.\n")) ; return (TRUE) ; } /* ========================================================================== */ /* === colamd_report ======================================================== */ /* ========================================================================== */ PUBLIC void colamd_report ( int stats [COLAMD_STATS] ) { print_report ("colamd", stats) ; } /* ========================================================================== */ /* === symamd_report ======================================================== */ /* ========================================================================== */ PUBLIC void symamd_report ( int stats [COLAMD_STATS] ) { print_report ("symamd", stats) ; } /* ========================================================================== */ /* === NON-USER-CALLABLE ROUTINES: ========================================== */ /* ========================================================================== */ /* There are no user-callable routines beyond this point in the file */ /* ========================================================================== */ /* === init_rows_cols ======================================================= */ /* ========================================================================== */ /* Takes the column form of the matrix in A and creates the row form of the matrix. Also, row and column attributes are stored in the Col and Row structs. If the columns are un-sorted or contain duplicate row indices, this routine will also sort and remove duplicate row indices from the column form of the matrix. Returns FALSE if the matrix is invalid, TRUE otherwise. Not user-callable. */ PRIVATE int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ Colamd_Row Row [], /* of size n_row+1 */ Colamd_Col Col [], /* of size n_col+1 */ int A [], /* row indices of A, of size Alen */ int p [], /* pointers to columns in A, of size n_col+1 */ int stats [COLAMD_STATS] /* colamd statistics */ ) { /* === Local variables ================================================== */ int col ; /* a column index */ int row ; /* a row index */ int *cp ; /* a column pointer */ int *cp_end ; /* a pointer to the end of a column */ int *rp ; /* a row pointer */ int *rp_end ; /* a pointer to the end of a row */ int last_row ; /* previous row */ /* === Initialize columns, and check column pointers ==================== */ for (col = 0 ; col < n_col ; col++) { Col [col].start = p [col] ; Col [col].length = p [col+1] - p [col] ; if (Col [col].length < 0) { /* column pointers must be non-decreasing */ stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; stats [COLAMD_INFO1] = col ; stats [COLAMD_INFO2] = Col [col].length ; DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; return (FALSE) ; } Col [col].shared1.thickness = 1 ; Col [col].shared2.score = 0 ; Col [col].shared3.prev = EMPTY ; Col [col].shared4.degree_next = EMPTY ; } /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ /* === Scan columns, compute row degrees, and check row indices ========= */ stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ for (row = 0 ; row < n_row ; row++) { Row [row].length = 0 ; Row [row].shared2.mark = -1 ; } for (col = 0 ; col < n_col ; col++) { last_row = -1 ; cp = &A [p [col]] ; cp_end = &A [p [col+1]] ; while (cp < cp_end) { row = *cp++ ; /* make sure row indices within range */ if (row < 0 || row >= n_row) { stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; stats [COLAMD_INFO1] = col ; stats [COLAMD_INFO2] = row ; stats [COLAMD_INFO3] = n_row ; DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; return (FALSE) ; } if (row <= last_row || Row [row].shared2.mark == col) { /* row index are unsorted or repeated (or both), thus col */ /* is jumbled. This is a notice, not an error condition. */ stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; stats [COLAMD_INFO1] = col ; stats [COLAMD_INFO2] = row ; (stats [COLAMD_INFO3]) ++ ; DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); } if (Row [row].shared2.mark != col) { Row [row].length++ ; } else { /* this is a repeated entry in the column, */ /* it will be removed */ Col [col].length-- ; } /* mark the row as having been seen in this column */ Row [row].shared2.mark = col ; last_row = row ; } } /* === Compute row pointers ============================================= */ /* row form of the matrix starts directly after the column */ /* form of matrix in A */ Row [0].start = p [n_col] ; Row [0].shared1.p = Row [0].start ; Row [0].shared2.mark = -1 ; for (row = 1 ; row < n_row ; row++) { Row [row].start = Row [row-1].start + Row [row-1].length ; Row [row].shared1.p = Row [row].start ; Row [row].shared2.mark = -1 ; } /* === Create row form ================================================== */ if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { /* if cols jumbled, watch for repeated row indices */ for (col = 0 ; col < n_col ; col++) { cp = &A [p [col]] ; cp_end = &A [p [col+1]] ; while (cp < cp_end) { row = *cp++ ; if (Row [row].shared2.mark != col) { A [(Row [row].shared1.p)++] = col ; Row [row].shared2.mark = col ; } } } } else { /* if cols not jumbled, we don't need the mark (this is faster) */ for (col = 0 ; col < n_col ; col++) { cp = &A [p [col]] ; cp_end = &A [p [col+1]] ; while (cp < cp_end) { A [(Row [*cp++].shared1.p)++] = col ; } } } /* === Clear the row marks and set row degrees ========================== */ for (row = 0 ; row < n_row ; row++) { Row [row].shared2.mark = 0 ; Row [row].shared1.degree = Row [row].length ; } /* === See if we need to re-create columns ============================== */ if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; #ifndef NDEBUG /* make sure column lengths are correct */ for (col = 0 ; col < n_col ; col++) { p [col] = Col [col].length ; } for (row = 0 ; row < n_row ; row++) { rp = &A [Row [row].start] ; rp_end = rp + Row [row].length ; while (rp < rp_end) { p [*rp++]-- ; } } for (col = 0 ; col < n_col ; col++) { ASSERT (p [col] == 0) ; } /* now p is all zero (different than when debugging is turned off) */ #endif /* NDEBUG */ /* === Compute col pointers ========================================= */ /* col form of the matrix starts at A [0]. */ /* Note, we may have a gap between the col form and the row */ /* form if there were duplicate entries, if so, it will be */ /* removed upon the first garbage collection */ Col [0].start = 0 ; p [0] = Col [0].start ; for (col = 1 ; col < n_col ; col++) { /* note that the lengths here are for pruned columns, i.e. */ /* no duplicate row indices will exist for these columns */ Col [col].start = Col [col-1].start + Col [col-1].length ; p [col] = Col [col].start ; } /* === Re-create col form =========================================== */ for (row = 0 ; row < n_row ; row++) { rp = &A [Row [row].start] ; rp_end = rp + Row [row].length ; while (rp < rp_end) { A [(p [*rp++])++] = row ; } } } /* === Done. Matrix is not (or no longer) jumbled ====================== */ return (TRUE) ; } /* ========================================================================== */ /* === init_scoring ========================================================= */ /* ========================================================================== */ /* Kills dense or empty columns and rows, calculates an initial score for each column, and places all columns in the degree lists. Not user-callable. */ PRIVATE void init_scoring ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ Colamd_Row Row [], /* of size n_row+1 */ Colamd_Col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ double knobs [COLAMD_KNOBS],/* parameters */ int *p_n_row2, /* number of non-dense, non-empty rows */ int *p_n_col2, /* number of non-dense, non-empty columns */ int *p_max_deg /* maximum row degree */ ) { /* === Local variables ================================================== */ int c ; /* a column index */ int r, row ; /* a row index */ int *cp ; /* a column pointer */ int deg ; /* degree of a row or column */ int *cp_end ; /* a pointer to the end of a column */ int *new_cp ; /* new column pointer */ int col_length ; /* length of pruned column */ int score ; /* current column score */ int n_col2 ; /* number of non-dense, non-empty columns */ int n_row2 ; /* number of non-dense, non-empty rows */ int dense_row_count ; /* remove rows with more entries than this */ int dense_col_count ; /* remove cols with more entries than this */ int min_score ; /* smallest column score */ int max_deg ; /* maximum row degree */ int next_col ; /* Used to add to degree list.*/ #ifndef NDEBUG int debug_count ; /* debug only. */ #endif /* NDEBUG */ /* === Extract knobs ==================================================== */ dense_row_count = MAX (0, MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ; dense_col_count = MAX (0, MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ; DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; max_deg = 0 ; n_col2 = n_col ; n_row2 = n_row ; /* === Kill empty columns =============================================== */ /* Put the empty columns at the end in their natural order, so that LU */ /* factorization can proceed as far as possible. */ for (c = n_col-1 ; c >= 0 ; c--) { deg = Col [c].length ; if (deg == 0) { /* this is a empty column, kill and order it last */ Col [c].shared2.order = --n_col2 ; KILL_PRINCIPAL_COL (c) ; } } DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense columns =============================================== */ /* Put the dense columns at the end, in their natural order */ for (c = n_col-1 ; c >= 0 ; c--) { /* skip any dead columns */ if (COL_IS_DEAD (c)) { continue ; } deg = Col [c].length ; if (deg > dense_col_count) { /* this is a dense column, kill and order it last */ Col [c].shared2.order = --n_col2 ; /* decrement the row degrees */ cp = &A [Col [c].start] ; cp_end = cp + Col [c].length ; while (cp < cp_end) { Row [*cp++].shared1.degree-- ; } KILL_PRINCIPAL_COL (c) ; } } DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense and empty rows ======================================== */ for (r = 0 ; r < n_row ; r++) { deg = Row [r].shared1.degree ; ASSERT (deg >= 0 && deg <= n_col) ; if (deg > dense_row_count || deg == 0) { /* kill a dense or empty row */ KILL_ROW (r) ; --n_row2 ; } else { /* keep track of max degree of remaining rows */ max_deg = MAX (max_deg, deg) ; } } DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; /* === Compute initial column scores ==================================== */ /* At this point the row degrees are accurate. They reflect the number */ /* of "live" (non-dense) columns in each row. No empty rows exist. */ /* Some "live" columns may contain only dead rows, however. These are */ /* pruned in the code below. */ /* now find the initial matlab score for each column */ for (c = n_col-1 ; c >= 0 ; c--) { /* skip dead column */ if (COL_IS_DEAD (c)) { continue ; } score = 0 ; cp = &A [Col [c].start] ; new_cp = cp ; cp_end = cp + Col [c].length ; while (cp < cp_end) { /* get a row */ row = *cp++ ; /* skip if dead */ if (ROW_IS_DEAD (row)) { continue ; } /* compact the column */ *new_cp++ = row ; /* add row's external degree */ score += Row [row].shared1.degree - 1 ; /* guard against integer overflow */ score = MIN (score, n_col) ; } /* determine pruned column length */ col_length = (int) (new_cp - &A [Col [c].start]) ; if (col_length == 0) { /* a newly-made null column (all rows in this col are "dense" */ /* and have already been killed) */ DEBUG2 (("Newly null killed: %d\n", c)) ; Col [c].shared2.order = --n_col2 ; KILL_PRINCIPAL_COL (c) ; } else { /* set column length and set score */ ASSERT (score >= 0) ; ASSERT (score <= n_col) ; Col [c].length = col_length ; Col [c].shared2.score = score ; } } DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", n_col-n_col2)) ; /* At this point, all empty rows and columns are dead. All live columns */ /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ /* yet). Rows may contain dead columns, but all live rows contain at */ /* least one live column. */ #ifndef NDEBUG debug_structures (n_row, n_col, Row, Col, A, n_col2) ; #endif /* NDEBUG */ /* === Initialize degree lists ========================================== */ #ifndef NDEBUG debug_count = 0 ; #endif /* NDEBUG */ /* clear the hash buckets */ for (c = 0 ; c <= n_col ; c++) { head [c] = EMPTY ; } min_score = n_col ; /* place in reverse order, so low column indices are at the front */ /* of the lists. This is to encourage natural tie-breaking */ for (c = n_col-1 ; c >= 0 ; c--) { /* only add principal columns to degree lists */ if (COL_IS_ALIVE (c)) { DEBUG4 (("place %d score %d minscore %d ncol %d\n", c, Col [c].shared2.score, min_score, n_col)) ; /* === Add columns score to DList =============================== */ score = Col [c].shared2.score ; ASSERT (min_score >= 0) ; ASSERT (min_score <= n_col) ; ASSERT (score >= 0) ; ASSERT (score <= n_col) ; ASSERT (head [score] >= EMPTY) ; /* now add this column to dList at proper score location */ next_col = head [score] ; Col [c].shared3.prev = EMPTY ; Col [c].shared4.degree_next = next_col ; /* if there already was a column with the same score, set its */ /* previous pointer to this new column */ if (next_col != EMPTY) { Col [next_col].shared3.prev = c ; } head [score] = c ; /* see if this score is less than current min */ min_score = MIN (min_score, score) ; #ifndef NDEBUG debug_count++ ; #endif /* NDEBUG */ } } #ifndef NDEBUG DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n", debug_count, n_col, n_col-debug_count)) ; ASSERT (debug_count == n_col2) ; debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; #endif /* NDEBUG */ /* === Return number of remaining columns, and max row degree =========== */ *p_n_col2 = n_col2 ; *p_n_row2 = n_row2 ; *p_max_deg = max_deg ; } /* ========================================================================== */ /* === find_ordering ======================================================== */ /* ========================================================================== */ /* Order the principal columns of the supercolumn form of the matrix (no supercolumns on input). Uses a minimum approximate column minimum degree ordering method. Not user-callable. */ PRIVATE int find_ordering /* return the number of garbage collections */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ int Alen, /* size of A, 2*nnz + n_col or larger */ Colamd_Row Row [], /* of size n_row+1 */ Colamd_Col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ int n_col2, /* Remaining columns to order */ int max_deg, /* Maximum row degree */ int pfree /* index of first free slot (2*nnz on entry) */ ) { /* === Local variables ================================================== */ int k ; /* current pivot ordering step */ int pivot_col ; /* current pivot column */ int *cp ; /* a column pointer */ int *rp ; /* a row pointer */ int pivot_row ; /* current pivot row */ int *new_cp ; /* modified column pointer */ int *new_rp ; /* modified row pointer */ int pivot_row_start ; /* pointer to start of pivot row */ int pivot_row_degree ; /* number of columns in pivot row */ int pivot_row_length ; /* number of supercolumns in pivot row */ int pivot_col_score ; /* score of pivot column */ int needed_memory ; /* free space needed for pivot row */ int *cp_end ; /* pointer to the end of a column */ int *rp_end ; /* pointer to the end of a row */ int row ; /* a row index */ int col ; /* a column index */ int max_score ; /* maximum possible score */ int cur_score ; /* score of current column */ unsigned int hash ; /* hash value for supernode detection */ int head_column ; /* head of hash bucket */ int first_col ; /* first column in hash bucket */ int tag_mark ; /* marker value for mark array */ int row_mark ; /* Row [row].shared2.mark */ int set_difference ; /* set difference size of row with pivot row */ int min_score ; /* smallest column score */ int col_thickness ; /* "thickness" (no. of columns in a supercol) */ int max_mark ; /* maximum value of tag_mark */ int pivot_col_thickness ; /* number of columns represented by pivot col */ int prev_col ; /* Used by Dlist operations. */ int next_col ; /* Used by Dlist operations. */ int ngarbage ; /* number of garbage collections performed */ #ifndef NDEBUG int debug_d ; /* debug loop counter */ int debug_step = 0 ; /* debug loop counter */ #endif /* NDEBUG */ /* === Initialization and clear mark ==================================== */ max_mark = INT_MAX - n_col ; /* INT_MAX defined in */ tag_mark = clear_mark (n_row, Row) ; min_score = 0 ; ngarbage = 0 ; DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; /* === Order the columns ================================================ */ for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) { #ifndef NDEBUG if (debug_step % 100 == 0) { DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; } else { DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; } debug_step++ ; debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2-k, max_deg) ; debug_matrix (n_row, n_col, Row, Col, A) ; #endif /* NDEBUG */ /* === Select pivot column, and order it ============================ */ /* make sure degree list isn't empty */ ASSERT (min_score >= 0) ; ASSERT (min_score <= n_col) ; ASSERT (head [min_score] >= EMPTY) ; #ifndef NDEBUG for (debug_d = 0 ; debug_d < min_score ; debug_d++) { ASSERT (head [debug_d] == EMPTY) ; } #endif /* NDEBUG */ /* get pivot column from head of minimum degree list */ while (head [min_score] == EMPTY && min_score < n_col) { min_score++ ; } pivot_col = head [min_score] ; ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; next_col = Col [pivot_col].shared4.degree_next ; head [min_score] = next_col ; if (next_col != EMPTY) { Col [next_col].shared3.prev = EMPTY ; } ASSERT (COL_IS_ALIVE (pivot_col)) ; DEBUG3 (("Pivot col: %d\n", pivot_col)) ; /* remember score for defrag check */ pivot_col_score = Col [pivot_col].shared2.score ; /* the pivot column is the kth column in the pivot order */ Col [pivot_col].shared2.order = k ; /* increment order count by column thickness */ pivot_col_thickness = Col [pivot_col].shared1.thickness ; k += pivot_col_thickness ; ASSERT (pivot_col_thickness > 0) ; /* === Garbage_collection, if necessary ============================= */ needed_memory = MIN (pivot_col_score, n_col - k) ; if (pfree + needed_memory >= Alen) { pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; ngarbage++ ; /* after garbage collection we will have enough */ ASSERT (pfree + needed_memory < Alen) ; /* garbage collection has wiped out the Row[].shared2.mark array */ tag_mark = clear_mark (n_row, Row) ; #ifndef NDEBUG debug_matrix (n_row, n_col, Row, Col, A) ; #endif /* NDEBUG */ } /* === Compute pivot row pattern ==================================== */ /* get starting location for this new merged row */ pivot_row_start = pfree ; /* initialize new row counts to zero */ pivot_row_degree = 0 ; /* tag pivot column as having been visited so it isn't included */ /* in merged pivot row */ Col [pivot_col].shared1.thickness = -pivot_col_thickness ; /* pivot row is the union of all rows in the pivot column pattern */ cp = &A [Col [pivot_col].start] ; cp_end = cp + Col [pivot_col].length ; while (cp < cp_end) { /* get a row */ row = *cp++ ; DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; /* skip if row is dead */ if (ROW_IS_DEAD (row)) { continue ; } rp = &A [Row [row].start] ; rp_end = rp + Row [row].length ; while (rp < rp_end) { /* get a column */ col = *rp++ ; /* add the column, if alive and untagged */ col_thickness = Col [col].shared1.thickness ; if (col_thickness > 0 && COL_IS_ALIVE (col)) { /* tag column in pivot row */ Col [col].shared1.thickness = -col_thickness ; ASSERT (pfree < Alen) ; /* place column in pivot row */ A [pfree++] = col ; pivot_row_degree += col_thickness ; } } } /* clear tag on pivot column */ Col [pivot_col].shared1.thickness = pivot_col_thickness ; max_deg = MAX (max_deg, pivot_row_degree) ; #ifndef NDEBUG DEBUG3 (("check2\n")) ; debug_mark (n_row, Row, tag_mark, max_mark) ; #endif /* NDEBUG */ /* === Kill all rows used to construct pivot row ==================== */ /* also kill pivot row, temporarily */ cp = &A [Col [pivot_col].start] ; cp_end = cp + Col [pivot_col].length ; while (cp < cp_end) { /* may be killing an already dead row */ row = *cp++ ; DEBUG3 (("Kill row in pivot col: %d\n", row)) ; KILL_ROW (row) ; } /* === Select a row index to use as the new pivot row =============== */ pivot_row_length = pfree - pivot_row_start ; if (pivot_row_length > 0) { /* pick the "pivot" row arbitrarily (first row in col) */ pivot_row = A [Col [pivot_col].start] ; DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; } else { /* there is no pivot row, since it is of zero length */ pivot_row = EMPTY ; ASSERT (pivot_row_length == 0) ; } ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; /* === Approximate degree computation =============================== */ /* Here begins the computation of the approximate degree. The column */ /* score is the sum of the pivot row "length", plus the size of the */ /* set differences of each row in the column minus the pattern of the */ /* pivot row itself. The column ("thickness") itself is also */ /* excluded from the column score (we thus use an approximate */ /* external degree). */ /* The time taken by the following code (compute set differences, and */ /* add them up) is proportional to the size of the data structure */ /* being scanned - that is, the sum of the sizes of each column in */ /* the pivot row. Thus, the amortized time to compute a column score */ /* is proportional to the size of that column (where size, in this */ /* context, is the column "length", or the number of row indices */ /* in that column). The number of row indices in a column is */ /* monotonically non-decreasing, from the length of the original */ /* column on input to colamd. */ /* === Compute set differences ====================================== */ DEBUG3 (("** Computing set differences phase. **\n")) ; /* pivot row is currently dead - it will be revived later. */ DEBUG3 (("Pivot row: ")) ; /* for each column in pivot row */ rp = &A [pivot_row_start] ; rp_end = rp + pivot_row_length ; while (rp < rp_end) { col = *rp++ ; ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; DEBUG3 (("Col: %d\n", col)) ; /* clear tags used to construct pivot row pattern */ col_thickness = -Col [col].shared1.thickness ; ASSERT (col_thickness > 0) ; Col [col].shared1.thickness = col_thickness ; /* === Remove column from degree list =========================== */ cur_score = Col [col].shared2.score ; prev_col = Col [col].shared3.prev ; next_col = Col [col].shared4.degree_next ; ASSERT (cur_score >= 0) ; ASSERT (cur_score <= n_col) ; ASSERT (cur_score >= EMPTY) ; if (prev_col == EMPTY) { head [cur_score] = next_col ; } else { Col [prev_col].shared4.degree_next = next_col ; } if (next_col != EMPTY) { Col [next_col].shared3.prev = prev_col ; } /* === Scan the column ========================================== */ cp = &A [Col [col].start] ; cp_end = cp + Col [col].length ; while (cp < cp_end) { /* get a row */ row = *cp++ ; row_mark = Row [row].shared2.mark ; /* skip if dead */ if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } ASSERT (row != pivot_row) ; set_difference = row_mark - tag_mark ; /* check if the row has been seen yet */ if (set_difference < 0) { ASSERT (Row [row].shared1.degree <= max_deg) ; set_difference = Row [row].shared1.degree ; } /* subtract column thickness from this row's set difference */ set_difference -= col_thickness ; ASSERT (set_difference >= 0) ; /* absorb this row if the set difference becomes zero */ if (set_difference == 0) { DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; KILL_ROW (row) ; } else { /* save the new mark */ Row [row].shared2.mark = set_difference + tag_mark ; } } } #ifndef NDEBUG debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2-k-pivot_row_degree, max_deg) ; #endif /* NDEBUG */ /* === Add up set differences for each column ======================= */ DEBUG3 (("** Adding set differences phase. **\n")) ; /* for each column in pivot row */ rp = &A [pivot_row_start] ; rp_end = rp + pivot_row_length ; while (rp < rp_end) { /* get a column */ col = *rp++ ; ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; hash = 0 ; cur_score = 0 ; cp = &A [Col [col].start] ; /* compact the column */ new_cp = cp ; cp_end = cp + Col [col].length ; DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; while (cp < cp_end) { /* get a row */ row = *cp++ ; ASSERT(row >= 0 && row < n_row) ; row_mark = Row [row].shared2.mark ; /* skip if dead */ if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } ASSERT (row_mark > tag_mark) ; /* compact the column */ *new_cp++ = row ; /* compute hash function */ hash += row ; /* add set difference */ cur_score += row_mark - tag_mark ; /* integer overflow... */ cur_score = MIN (cur_score, n_col) ; } /* recompute the column's length */ Col [col].length = (int) (new_cp - &A [Col [col].start]) ; /* === Further mass elimination ================================= */ if (Col [col].length == 0) { DEBUG4 (("further mass elimination. Col: %d\n", col)) ; /* nothing left but the pivot row in this column */ KILL_PRINCIPAL_COL (col) ; pivot_row_degree -= Col [col].shared1.thickness ; ASSERT (pivot_row_degree >= 0) ; /* order it */ Col [col].shared2.order = k ; /* increment order count by column thickness */ k += Col [col].shared1.thickness ; } else { /* === Prepare for supercolumn detection ==================== */ DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; /* save score so far */ Col [col].shared2.score = cur_score ; /* add column to hash table, for supercolumn detection */ hash %= n_col + 1 ; DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; ASSERT (hash <= n_col) ; head_column = head [hash] ; if (head_column > EMPTY) { /* degree list "hash" is non-empty, use prev (shared3) of */ /* first column in degree list as head of hash bucket */ first_col = Col [head_column].shared3.headhash ; Col [head_column].shared3.headhash = col ; } else { /* degree list "hash" is empty, use head as hash bucket */ first_col = - (head_column + 2) ; head [hash] = - (col + 2) ; } Col [col].shared4.hash_next = first_col ; /* save hash function in Col [col].shared3.hash */ Col [col].shared3.hash = (int) hash ; ASSERT (COL_IS_ALIVE (col)) ; } } /* The approximate external column degree is now computed. */ /* === Supercolumn detection ======================================== */ DEBUG3 (("** Supercolumn detection phase. **\n")) ; detect_super_cols ( #ifndef NDEBUG n_col, Row, #endif /* NDEBUG */ Col, A, head, pivot_row_start, pivot_row_length) ; /* === Kill the pivotal column ====================================== */ KILL_PRINCIPAL_COL (pivot_col) ; /* === Clear mark =================================================== */ tag_mark += (max_deg + 1) ; if (tag_mark >= max_mark) { DEBUG2 (("clearing tag_mark\n")) ; tag_mark = clear_mark (n_row, Row) ; } #ifndef NDEBUG DEBUG3 (("check3\n")) ; debug_mark (n_row, Row, tag_mark, max_mark) ; #endif /* NDEBUG */ /* === Finalize the new pivot row, and column scores ================ */ DEBUG3 (("** Finalize scores phase. **\n")) ; /* for each column in pivot row */ rp = &A [pivot_row_start] ; /* compact the pivot row */ new_rp = rp ; rp_end = rp + pivot_row_length ; while (rp < rp_end) { col = *rp++ ; /* skip dead columns */ if (COL_IS_DEAD (col)) { continue ; } *new_rp++ = col ; /* add new pivot row to column */ A [Col [col].start + (Col [col].length++)] = pivot_row ; /* retrieve score so far and add on pivot row's degree. */ /* (we wait until here for this in case the pivot */ /* row's degree was reduced due to mass elimination). */ cur_score = Col [col].shared2.score + pivot_row_degree ; /* calculate the max possible score as the number of */ /* external columns minus the 'k' value minus the */ /* columns thickness */ max_score = n_col - k - Col [col].shared1.thickness ; /* make the score the external degree of the union-of-rows */ cur_score -= Col [col].shared1.thickness ; /* make sure score is less or equal than the max score */ cur_score = MIN (cur_score, max_score) ; ASSERT (cur_score >= 0) ; /* store updated score */ Col [col].shared2.score = cur_score ; /* === Place column back in degree list ========================= */ ASSERT (min_score >= 0) ; ASSERT (min_score <= n_col) ; ASSERT (cur_score >= 0) ; ASSERT (cur_score <= n_col) ; ASSERT (head [cur_score] >= EMPTY) ; next_col = head [cur_score] ; Col [col].shared4.degree_next = next_col ; Col [col].shared3.prev = EMPTY ; if (next_col != EMPTY) { Col [next_col].shared3.prev = col ; } head [cur_score] = col ; /* see if this score is less than current min */ min_score = MIN (min_score, cur_score) ; } #ifndef NDEBUG debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2-k, max_deg) ; #endif /* NDEBUG */ /* === Resurrect the new pivot row ================================== */ if (pivot_row_degree > 0) { /* update pivot row length to reflect any cols that were killed */ /* during super-col detection and mass elimination */ Row [pivot_row].start = pivot_row_start ; Row [pivot_row].length = (int) (new_rp - &A[pivot_row_start]) ; Row [pivot_row].shared1.degree = pivot_row_degree ; Row [pivot_row].shared2.mark = 0 ; /* pivot row is no longer dead */ } } /* === All principal columns have now been ordered ====================== */ return (ngarbage) ; } /* ========================================================================== */ /* === order_children ======================================================= */ /* ========================================================================== */ /* The find_ordering routine has ordered all of the principal columns (the representatives of the supercolumns). The non-principal columns have not yet been ordered. This routine orders those columns by walking up the parent tree (a column is a child of the column which absorbed it). The final permutation vector is then placed in p [0 ... n_col-1], with p [0] being the first column, and p [n_col-1] being the last. It doesn't look like it at first glance, but be assured that this routine takes time linear in the number of columns. Although not immediately obvious, the time taken by this routine is O (n_col), that is, linear in the number of columns. Not user-callable. */ PRIVATE void order_children ( /* === Parameters ======================================================= */ int n_col, /* number of columns of A */ Colamd_Col Col [], /* of size n_col+1 */ int p [] /* p [0 ... n_col-1] is the column permutation*/ ) { /* === Local variables ================================================== */ int i ; /* loop counter for all columns */ int c ; /* column index */ int parent ; /* index of column's parent */ int order ; /* column's order */ /* === Order each non-principal column ================================== */ for (i = 0 ; i < n_col ; i++) { /* find an un-ordered non-principal column */ ASSERT (COL_IS_DEAD (i)) ; if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY) { parent = i ; /* once found, find its principal parent */ do { parent = Col [parent].shared1.parent ; } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; /* now, order all un-ordered non-principal columns along path */ /* to this parent. collapse tree at the same time */ c = i ; /* get order of parent */ order = Col [parent].shared2.order ; do { ASSERT (Col [c].shared2.order == EMPTY) ; /* order this column */ Col [c].shared2.order = order++ ; /* collaps tree */ Col [c].shared1.parent = parent ; /* get immediate parent of this column */ c = Col [c].shared1.parent ; /* continue until we hit an ordered column. There are */ /* guarranteed not to be anymore unordered columns */ /* above an ordered column */ } while (Col [c].shared2.order == EMPTY) ; /* re-order the super_col parent to largest order for this group */ Col [parent].shared2.order = order ; } } /* === Generate the permutation ========================================= */ for (c = 0 ; c < n_col ; c++) { p [Col [c].shared2.order] = c ; } } /* ========================================================================== */ /* === detect_super_cols ==================================================== */ /* ========================================================================== */ /* Detects supercolumns by finding matches between columns in the hash buckets. Check amongst columns in the set A [row_start ... row_start + row_length-1]. The columns under consideration are currently *not* in the degree lists, and have already been placed in the hash buckets. The hash bucket for columns whose hash function is equal to h is stored as follows: if head [h] is >= 0, then head [h] contains a degree list, so: head [h] is the first column in degree bucket h. Col [head [h]].headhash gives the first column in hash bucket h. otherwise, the degree list is empty, and: -(head [h] + 2) is the first column in hash bucket h. For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous column" pointer. Col [c].shared3.hash is used instead as the hash number for that column. The value of Col [c].shared4.hash_next is the next column in the same hash bucket. Assuming no, or "few" hash collisions, the time taken by this routine is linear in the sum of the sizes (lengths) of each column whose score has just been computed in the approximate degree computation. Not user-callable. */ PRIVATE void detect_super_cols ( /* === Parameters ======================================================= */ #ifndef NDEBUG /* these two parameters are only needed when debugging is enabled: */ int n_col, /* number of columns of A */ Colamd_Row Row [], /* of size n_row+1 */ #endif /* NDEBUG */ Colamd_Col Col [], /* of size n_col+1 */ int A [], /* row indices of A */ int head [], /* head of degree lists and hash buckets */ int row_start, /* pointer to set of columns to check */ int row_length /* number of columns to check */ ) { /* === Local variables ================================================== */ int hash ; /* hash value for a column */ int *rp ; /* pointer to a row */ int c ; /* a column index */ int super_c ; /* column index of the column to absorb into */ int *cp1 ; /* column pointer for column super_c */ int *cp2 ; /* column pointer for column c */ int length ; /* length of column super_c */ int prev_c ; /* column preceding c in hash bucket */ int i ; /* loop counter */ int *rp_end ; /* pointer to the end of the row */ int col ; /* a column index in the row to check */ int head_column ; /* first column in hash bucket or degree list */ int first_col ; /* first column in hash bucket */ /* === Consider each column in the row ================================== */ rp = &A [row_start] ; rp_end = rp + row_length ; while (rp < rp_end) { col = *rp++ ; if (COL_IS_DEAD (col)) { continue ; } /* get hash number for this column */ hash = Col [col].shared3.hash ; ASSERT (hash <= n_col) ; /* === Get the first column in this hash bucket ===================== */ head_column = head [hash] ; if (head_column > EMPTY) { first_col = Col [head_column].shared3.headhash ; } else { first_col = - (head_column + 2) ; } /* === Consider each column in the hash bucket ====================== */ for (super_c = first_col ; super_c != EMPTY ; super_c = Col [super_c].shared4.hash_next) { ASSERT (COL_IS_ALIVE (super_c)) ; ASSERT (Col [super_c].shared3.hash == hash) ; length = Col [super_c].length ; /* prev_c is the column preceding column c in the hash bucket */ prev_c = super_c ; /* === Compare super_c with all columns after it ================ */ for (c = Col [super_c].shared4.hash_next ; c != EMPTY ; c = Col [c].shared4.hash_next) { ASSERT (c != super_c) ; ASSERT (COL_IS_ALIVE (c)) ; ASSERT (Col [c].shared3.hash == hash) ; /* not identical if lengths or scores are different */ if (Col [c].length != length || Col [c].shared2.score != Col [super_c].shared2.score) { prev_c = c ; continue ; } /* compare the two columns */ cp1 = &A [Col [super_c].start] ; cp2 = &A [Col [c].start] ; for (i = 0 ; i < length ; i++) { /* the columns are "clean" (no dead rows) */ ASSERT (ROW_IS_ALIVE (*cp1)) ; ASSERT (ROW_IS_ALIVE (*cp2)) ; /* row indices will same order for both supercols, */ /* no gather scatter nessasary */ if (*cp1++ != *cp2++) { break ; } } /* the two columns are different if the for-loop "broke" */ if (i != length) { prev_c = c ; continue ; } /* === Got it! two columns are identical =================== */ ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; Col [super_c].shared1.thickness += Col [c].shared1.thickness ; Col [c].shared1.parent = super_c ; KILL_NON_PRINCIPAL_COL (c) ; /* order c later, in order_children() */ Col [c].shared2.order = EMPTY ; /* remove c from hash bucket */ Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; } } /* === Empty this hash bucket ======================================= */ if (head_column > EMPTY) { /* corresponding degree list "hash" is not empty */ Col [head_column].shared3.headhash = EMPTY ; } else { /* corresponding degree list "hash" is empty */ head [hash] = EMPTY ; } } } /* ========================================================================== */ /* === garbage_collection =================================================== */ /* ========================================================================== */ /* Defragments and compacts columns and rows in the workspace A. Used when all avaliable memory has been used while performing row merging. Returns the index of the first free position in A, after garbage collection. The time taken by this routine is linear is the size of the array A, which is itself linear in the number of nonzeros in the input matrix. Not user-callable. */ PRIVATE int garbage_collection /* returns the new value of pfree */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows */ int n_col, /* number of columns */ Colamd_Row Row [], /* row info */ Colamd_Col Col [], /* column info */ int A [], /* A [0 ... Alen-1] holds the matrix */ int *pfree /* &A [0] ... pfree is in use */ ) { /* === Local variables ================================================== */ int *psrc ; /* source pointer */ int *pdest ; /* destination pointer */ int j ; /* counter */ int r ; /* a row index */ int c ; /* a column index */ int length ; /* length of a row or column */ #ifndef NDEBUG int debug_rows ; DEBUG2 (("Defrag..\n")) ; for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ; debug_rows = 0 ; #endif /* NDEBUG */ /* === Defragment the columns =========================================== */ pdest = &A[0] ; for (c = 0 ; c < n_col ; c++) { if (COL_IS_ALIVE (c)) { psrc = &A [Col [c].start] ; /* move and compact the column */ ASSERT (pdest <= psrc) ; Col [c].start = (int) (pdest - &A [0]) ; length = Col [c].length ; for (j = 0 ; j < length ; j++) { r = *psrc++ ; if (ROW_IS_ALIVE (r)) { *pdest++ = r ; } } Col [c].length = (int) (pdest - &A [Col [c].start]) ; } } /* === Prepare to defragment the rows =================================== */ for (r = 0 ; r < n_row ; r++) { if (ROW_IS_ALIVE (r)) { if (Row [r].length == 0) { /* this row is of zero length. cannot compact it, so kill it */ DEBUG3 (("Defrag row kill\n")) ; KILL_ROW (r) ; } else { /* save first column index in Row [r].shared2.first_column */ psrc = &A [Row [r].start] ; Row [r].shared2.first_column = *psrc ; ASSERT (ROW_IS_ALIVE (r)) ; /* flag the start of the row with the one's complement of row */ *psrc = ONES_COMPLEMENT (r) ; #ifndef NDEBUG debug_rows++ ; #endif /* NDEBUG */ } } } /* === Defragment the rows ============================================== */ psrc = pdest ; while (psrc < pfree) { /* find a negative number ... the start of a row */ if (*psrc++ < 0) { psrc-- ; /* get the row index */ r = ONES_COMPLEMENT (*psrc) ; ASSERT (r >= 0 && r < n_row) ; /* restore first column index */ *psrc = Row [r].shared2.first_column ; ASSERT (ROW_IS_ALIVE (r)) ; /* move and compact the row */ ASSERT (pdest <= psrc) ; Row [r].start = (int) (pdest - &A [0]) ; length = Row [r].length ; for (j = 0 ; j < length ; j++) { c = *psrc++ ; if (COL_IS_ALIVE (c)) { *pdest++ = c ; } } Row [r].length = (int) (pdest - &A [Row [r].start]) ; #ifndef NDEBUG debug_rows-- ; #endif /* NDEBUG */ } } /* ensure we found all the rows */ ASSERT (debug_rows == 0) ; /* === Return the new value of pfree ==================================== */ return ((int) (pdest - &A [0])) ; } /* ========================================================================== */ /* === clear_mark =========================================================== */ /* ========================================================================== */ /* Clears the Row [].shared2.mark array, and returns the new tag_mark. Return value is the new tag_mark. Not user-callable. */ PRIVATE int clear_mark /* return the new value for tag_mark */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows in A */ Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ ) { /* === Local variables ================================================== */ int r ; for (r = 0 ; r < n_row ; r++) { if (ROW_IS_ALIVE (r)) { Row [r].shared2.mark = 0 ; } } return (1) ; } /* ========================================================================== */ /* === print_report ========================================================= */ /* ========================================================================== */ PRIVATE void print_report ( char *method, int stats [COLAMD_STATS] ) { int i1, i2, i3 ; if (!stats) { PRINTF ("%s: No statistics available.\n", method) ; return ; } i1 = stats [COLAMD_INFO1] ; i2 = stats [COLAMD_INFO2] ; i3 = stats [COLAMD_INFO3] ; if (stats [COLAMD_STATUS] >= 0) { PRINTF ("%s: OK. ", method) ; } else { PRINTF ("%s: ERROR. ", method) ; } switch (stats [COLAMD_STATUS]) { case COLAMD_OK_BUT_JUMBLED: PRINTF ("Matrix has unsorted or duplicate row indices.\n") ; PRINTF ("%s: number of duplicate or out-of-order row indices: %d\n", method, i3) ; PRINTF ("%s: last seen duplicate or out-of-order row index: %d\n", method, INDEX (i2)) ; PRINTF ("%s: last seen in column: %d", method, INDEX (i1)) ; /* no break - fall through to next case instead */ case COLAMD_OK: PRINTF ("\n") ; PRINTF ("%s: number of dense or empty rows ignored: %d\n", method, stats [COLAMD_DENSE_ROW]) ; PRINTF ("%s: number of dense or empty columns ignored: %d\n", method, stats [COLAMD_DENSE_COL]) ; PRINTF ("%s: number of garbage collections performed: %d\n", method, stats [COLAMD_DEFRAG_COUNT]) ; break ; case COLAMD_ERROR_A_not_present: PRINTF ("Array A (row indices of matrix) not present.\n") ; break ; case COLAMD_ERROR_p_not_present: PRINTF ("Array p (column pointers for matrix) not present.\n") ; break ; case COLAMD_ERROR_nrow_negative: PRINTF ("Invalid number of rows (%d).\n", i1) ; break ; case COLAMD_ERROR_ncol_negative: PRINTF ("Invalid number of columns (%d).\n", i1) ; break ; case COLAMD_ERROR_nnz_negative: PRINTF ("Invalid number of nonzero entries (%d).\n", i1) ; break ; case COLAMD_ERROR_p0_nonzero: PRINTF ("Invalid column pointer, p [0] = %d, must be zero.\n", i1) ; break ; case COLAMD_ERROR_A_too_small: PRINTF ("Array A too small.\n") ; PRINTF (" Need Alen >= %d, but given only Alen = %d.\n", i1, i2) ; break ; case COLAMD_ERROR_col_length_negative: PRINTF ("Column %d has a negative number of nonzero entries (%d).\n", INDEX (i1), i2) ; break ; case COLAMD_ERROR_row_index_out_of_bounds: PRINTF ("Row index (row %d) out of bounds (%d to %d) in column %d.\n", INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1)) ; break ; case COLAMD_ERROR_out_of_memory: PRINTF ("Out of memory.\n") ; break ; case COLAMD_ERROR_internal_error: /* if this happens, there is a bug in the code */ PRINTF ("Internal error! Please contact authors (davis@cise.ufl.edu).\n") ; break ; } } /* ========================================================================== */ /* === colamd debugging routines ============================================ */ /* ========================================================================== */ /* When debugging is disabled, the remainder of this file is ignored. */ #ifndef NDEBUG /* ========================================================================== */ /* === debug_structures ===================================================== */ /* ========================================================================== */ /* At this point, all empty rows and columns are dead. All live columns are "clean" (containing no dead rows) and simplicial (no supercolumns yet). Rows may contain dead columns, but all live rows contain at least one live column. */ PRIVATE void debug_structures ( /* === Parameters ======================================================= */ int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [], int n_col2 ) { /* === Local variables ================================================== */ int i ; int c ; int *cp ; int *cp_end ; int len ; int score ; int r ; int *rp ; int *rp_end ; int deg ; /* === Check A, Row, and Col ============================================ */ for (c = 0 ; c < n_col ; c++) { if (COL_IS_ALIVE (c)) { len = Col [c].length ; score = Col [c].shared2.score ; DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; ASSERT (len > 0) ; ASSERT (score >= 0) ; ASSERT (Col [c].shared1.thickness == 1) ; cp = &A [Col [c].start] ; cp_end = cp + len ; while (cp < cp_end) { r = *cp++ ; ASSERT (ROW_IS_ALIVE (r)) ; } } else { i = Col [c].shared2.order ; ASSERT (i >= n_col2 && i < n_col) ; } } for (r = 0 ; r < n_row ; r++) { if (ROW_IS_ALIVE (r)) { i = 0 ; len = Row [r].length ; deg = Row [r].shared1.degree ; ASSERT (len > 0) ; ASSERT (deg > 0) ; rp = &A [Row [r].start] ; rp_end = rp + len ; while (rp < rp_end) { c = *rp++ ; if (COL_IS_ALIVE (c)) { i++ ; } } ASSERT (i > 0) ; } } } /* ========================================================================== */ /* === debug_deg_lists ====================================================== */ /* ========================================================================== */ /* Prints the contents of the degree lists. Counts the number of columns in the degree list and compares it to the total it should have. Also checks the row degrees. */ PRIVATE void debug_deg_lists ( /* === Parameters ======================================================= */ int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int head [], int min_score, int should, int max_deg ) { /* === Local variables ================================================== */ int deg ; int col ; int have ; int row ; /* === Check the degree lists =========================================== */ if (n_col > 10000 && colamd_debug <= 0) { return ; } have = 0 ; DEBUG4 (("Degree lists: %d\n", min_score)) ; for (deg = 0 ; deg <= n_col ; deg++) { col = head [deg] ; if (col == EMPTY) { continue ; } DEBUG4 (("%d:", deg)) ; while (col != EMPTY) { DEBUG4 ((" %d", col)) ; have += Col [col].shared1.thickness ; ASSERT (COL_IS_ALIVE (col)) ; col = Col [col].shared4.degree_next ; } DEBUG4 (("\n")) ; } DEBUG4 (("should %d have %d\n", should, have)) ; ASSERT (should == have) ; /* === Check the row degrees ============================================ */ if (n_row > 10000 && colamd_debug <= 0) { return ; } for (row = 0 ; row < n_row ; row++) { if (ROW_IS_ALIVE (row)) { ASSERT (Row [row].shared1.degree <= max_deg) ; } } } /* ========================================================================== */ /* === debug_mark =========================================================== */ /* ========================================================================== */ /* Ensures that the tag_mark is less that the maximum and also ensures that each entry in the mark array is less than the tag mark. */ PRIVATE void debug_mark ( /* === Parameters ======================================================= */ int n_row, Colamd_Row Row [], int tag_mark, int max_mark ) { /* === Local variables ================================================== */ int r ; /* === Check the Row marks ============================================== */ ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; if (n_row > 10000 && colamd_debug <= 0) { return ; } for (r = 0 ; r < n_row ; r++) { ASSERT (Row [r].shared2.mark < tag_mark) ; } } /* ========================================================================== */ /* === debug_matrix ========================================================= */ /* ========================================================================== */ /* Prints out the contents of the columns and the rows. */ PRIVATE void debug_matrix ( /* === Parameters ======================================================= */ int n_row, int n_col, Colamd_Row Row [], Colamd_Col Col [], int A [] ) { /* === Local variables ================================================== */ int r ; int c ; int *rp ; int *rp_end ; int *cp ; int *cp_end ; /* === Dump the rows and columns of the matrix ========================== */ if (colamd_debug < 3) { return ; } DEBUG3 (("DUMP MATRIX:\n")) ; for (r = 0 ; r < n_row ; r++) { DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ; if (ROW_IS_DEAD (r)) { continue ; } DEBUG3 (("start %d length %d degree %d\n", Row [r].start, Row [r].length, Row [r].shared1.degree)) ; rp = &A [Row [r].start] ; rp_end = rp + Row [r].length ; while (rp < rp_end) { c = *rp++ ; DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ; } } for (c = 0 ; c < n_col ; c++) { DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ; if (COL_IS_DEAD (c)) { continue ; } DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", Col [c].start, Col [c].length, Col [c].shared1.thickness, Col [c].shared2.score)) ; cp = &A [Col [c].start] ; cp_end = cp + Col [c].length ; while (cp < cp_end) { r = *cp++ ; DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ; } } } PRIVATE void colamd_get_debug ( char *method ) { colamd_debug = 0 ; /* no debug printing */ /* get "D" environment variable, which gives the debug printing level */ if (getenv ("D")) { colamd_debug = atoi (getenv ("D")) ; } DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", method, colamd_debug)) ; } #endif /* NDEBUG */ SHAR_EOF fi # end of overwriting check if test -f 'colamd.h' then echo shar: will not over-write existing file "'colamd.h'" else cat << "SHAR_EOF" > 'colamd.h' /* ========================================================================== */ /* === colamd/symamd prototypes and definitions ============================= */ /* ========================================================================== */ /* You must include this file (colamd.h) in any routine that uses colamd, symamd, or the related macros and definitions. Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Notice: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. Permission is hereby granted to use, copy, modify, and/or distribute this program, provided that the Copyright, this License, and the Availability of the original version is retained on all copies and made accessible to the end-user of any code or package that includes COLAMD or any modified version of COLAMD. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h file. It is required by the colamd.c, colamdmex.c, and symamdmex.c files, and by any C code that calls the routines whose prototypes are listed below, or that uses the colamd/symamd definitions listed below. */ #ifndef COLAMD_H #define COLAMD_H /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include /* ========================================================================== */ /* === Knob and statistics definitions ====================================== */ /* ========================================================================== */ /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ #define COLAMD_KNOBS 20 /* number of output statistics. Only stats [0..6] are currently used. */ #define COLAMD_STATS 20 /* knobs [0] and stats [0]: dense row knob and output statistic. */ #define COLAMD_DENSE_ROW 0 /* knobs [1] and stats [1]: dense column knob and output statistic. */ #define COLAMD_DENSE_COL 1 /* stats [2]: memory defragmentation count output statistic */ #define COLAMD_DEFRAG_COUNT 2 /* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ #define COLAMD_STATUS 3 /* stats [4..6]: error info, or info on jumbled columns */ #define COLAMD_INFO1 4 #define COLAMD_INFO2 5 #define COLAMD_INFO3 6 /* error codes returned in stats [3]: */ #define COLAMD_OK (0) #define COLAMD_OK_BUT_JUMBLED (1) #define COLAMD_ERROR_A_not_present (-1) #define COLAMD_ERROR_p_not_present (-2) #define COLAMD_ERROR_nrow_negative (-3) #define COLAMD_ERROR_ncol_negative (-4) #define COLAMD_ERROR_nnz_negative (-5) #define COLAMD_ERROR_p0_nonzero (-6) #define COLAMD_ERROR_A_too_small (-7) #define COLAMD_ERROR_col_length_negative (-8) #define COLAMD_ERROR_row_index_out_of_bounds (-9) #define COLAMD_ERROR_out_of_memory (-10) #define COLAMD_ERROR_internal_error (-999) /* ========================================================================== */ /* === Row and Column structures ============================================ */ /* ========================================================================== */ /* User code that makes use of the colamd/symamd routines need not directly */ /* reference these structures. They are used only for the COLAMD_RECOMMENDED */ /* macro. */ typedef struct Colamd_Col_struct { int start ; /* index for A of first row in this column, or DEAD */ /* if column is dead */ int length ; /* number of rows in this column */ union { int thickness ; /* number of original columns represented by this */ /* col, if the column is alive */ int parent ; /* parent in parent tree super-column structure, if */ /* the column is dead */ } shared1 ; union { int score ; /* the score used to maintain heap, if col is alive */ int order ; /* pivot ordering of this column, if col is dead */ } shared2 ; union { int headhash ; /* head of a hash bucket, if col is at the head of */ /* a degree list */ int hash ; /* hash value, if col is not in a degree list */ int prev ; /* previous column in degree list, if col is in a */ /* degree list (but not at the head of a degree list) */ } shared3 ; union { int degree_next ; /* next column, if col is in a degree list */ int hash_next ; /* next column, if col is in a hash list */ } shared4 ; } Colamd_Col ; typedef struct Colamd_Row_struct { int start ; /* index for A of first col in this row */ int length ; /* number of principal columns in this row */ union { int degree ; /* number of principal & non-principal columns in row */ int p ; /* used as a row pointer in init_rows_cols () */ } shared1 ; union { int mark ; /* for computing set differences and marking dead rows*/ int first_column ;/* first column in row (used in garbage collection) */ } shared2 ; } Colamd_Row ; /* ========================================================================== */ /* === Colamd recommended memory size ======================================= */ /* ========================================================================== */ /* The recommended length Alen of the array A passed to colamd is given by the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any argument is negative. 2*nnz space is required for the row and column indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is required for the Col and Row arrays, respectively, which are internal to colamd. An additional n_col space is the minimal amount of "elbow room", and nnz/5 more space is recommended for run time efficiency. This macro is not needed when using symamd. Explicit typecast to int added Sept. 23, 2002, COLAMD version 2.2, to avoid gcc -pedantic warning messages. */ #define COLAMD_C(n_col) ((int) (((n_col) + 1) * sizeof (Colamd_Col) / sizeof (int))) #define COLAMD_R(n_row) ((int) (((n_row) + 1) * sizeof (Colamd_Row) / sizeof (int))) #define COLAMD_RECOMMENDED(nnz, n_row, n_col) \ ( \ ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) \ ? \ (-1) \ : \ (2 * (nnz) + COLAMD_C (n_col) + COLAMD_R (n_row) + (n_col) + ((nnz) / 5)) \ ) /* ========================================================================== */ /* === Prototypes of user-callable routines ================================= */ /* ========================================================================== */ int colamd_recommended /* returns recommended value of Alen, */ /* or (-1) if input arguments are erroneous */ ( int nnz, /* nonzeros in A */ int n_row, /* number of rows in A */ int n_col /* number of columns in A */ ) ; void colamd_set_defaults /* sets default parameters */ ( /* knobs argument is modified on output */ double knobs [COLAMD_KNOBS] /* parameter settings for colamd */ ) ; int colamd /* returns (1) if successful, (0) otherwise*/ ( /* A and p arguments are modified on output */ int n_row, /* number of rows in A */ int n_col, /* number of columns in A */ int Alen, /* size of the array A */ int A [], /* row indices of A, of size Alen */ int p [], /* column pointers of A, of size n_col+1 */ double knobs [COLAMD_KNOBS],/* parameter settings for colamd */ int stats [COLAMD_STATS] /* colamd output statistics and error codes */ ) ; int symamd /* return (1) if OK, (0) otherwise */ ( int n, /* number of rows and columns of A */ int A [], /* row indices of A */ int p [], /* column pointers of A */ int perm [], /* output permutation, size n_col+1 */ double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */ int stats [COLAMD_STATS], /* output statistics and error codes */ void * (*allocate) (size_t, size_t), /* pointer to calloc (ANSI C) or */ /* mxCalloc (for MATLAB mexFunction) */ void (*release) (void *) /* pointer to free (ANSI C) or */ /* mxFree (for MATLAB mexFunction) */ ) ; void colamd_report ( int stats [COLAMD_STATS] ) ; void symamd_report ( int stats [COLAMD_STATS] ) ; #endif /* COLAMD_H */ SHAR_EOF fi # end of overwriting check cd .. cd .. if test -f 'ChangeLog' then echo shar: will not over-write existing file "'ChangeLog'" else cat << "SHAR_EOF" > 'ChangeLog' Changes from Version 2.2 to 2.3 (Sept. 8, 2003) * removed the call to the MATLAB spparms ('spumoni') function. This can take a lot of time if you are ordering many small matrices. Only affects the MATLAB interface (colamdmex.c, symamdmex.c, colamdtestmex.c, and symamdtestmex.c). The usage of the optional 2nd argument to the colamd and symamd mexFunctions was changed accordingly. Changes from Version 2.1 to 2.2 (Sept. 23, 2002) * extensive testing routines added (colamd_test.m, colamdtestmex.c, and symamdtestmex.c), and the Makefile modified accordingly. * a few typos in the comments corrected * use of the MATLAB "flops" command removed from colamd_demo, and an m-file routine luflops.m added. * an explicit typecast from unsigned to int added, for COLAMD_C and COLAMD_R in colamd.h. * #include added to colamd_example.c Changes from Version 2.0 to 2.1 (May 4, 2001) * TRUE and FALSE are predefined on some systems, so they are defined here only if not already defined. * web site changed * UNIX Makefile modified, to handle the case if "." is not in your path. Changes from Version 1.0 to 2.0 (January 31, 2000) No bugs were found in version 1.1. These changes merely add new functionality. * added the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. * moved the output statistics, from A, to a separate output argument. The arguments changed for the C-callable routines. * added colamd_report and symamd_report. * added a C-callable symamd routine. Formerly, symamd was only available as a mexFunction from MATLAB. * added error-checking to symamd. Formerly, it assumed its input was error-free. * added the optional stats and knobs arguments to the symamd mexFunction * deleted colamd_help. A help message is still available from "help colamd" and "help symamd" in MATLAB. * deleted colamdtree.m and symamdtree.m. Now, colamd.m and symamd.m also do the elimination tree post-ordering. The Version 1.1 colamd and symamd mexFunctions, which do not do the post- ordering, are now visible as colamdmex and symamdmex from MATLAB. Essentialy, the post-ordering is now the default behavior of colamd.m and symamd.m, to match the behavior of colmmd and symmmd. The post-ordering is only available in the MATLAB interface, not the C-callable interface. * made a slight change to the dense row/column detection in symamd, to match the stated specifications. SHAR_EOF fi # end of overwriting check if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << "SHAR_EOF" > 'Makefile' colamd_example: colamd_example.c colamd.c colamd.h cc -O -o colamd_example colamd_example.c colamd.c - ./colamd_example clean: - rm *.o colamd_example - rm colamdmex.mex* symamdmex.mex* - rm colamdtestmex.mex* symamdtestmex.mex* # Compiles the MATLAB-callable routines matlab: colamdmex.c colamd.c colamd.h mex -O colamdmex.c colamd.c mex -O symamdmex.c colamd.c # Compiles the extensive test code test: matlab colamdmex.c colamd.c colamd.h mex -O colamdtestmex.c colamd.c mex -O symamdtestmex.c colamd.c SHAR_EOF fi # end of overwriting check if test ! -d 'Matlab' then mkdir 'Matlab' fi cd 'Matlab' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'colamd_demo.m' then echo shar: will not over-write existing file "'colamd_demo.m'" else cat << "SHAR_EOF" > 'colamd_demo.m' % Demo for colamd: column approximate minimum degree ordering algorithm. % % The following m-files and mexFunctions provide alternative sparse matrix % ordering methods for MATLAB. They are typically faster (sometimes much % faster) and typically provide better orderings than their MATLAB counterparts: % % colamd a replacement for colmmd. % % Typical usage: p = colamd (A) ; % % symamd a replacement for symmmd. Based on colamd. % % Typical usage: p = symamd (A) ; % % For a description of the methods used, see the colamd.c file. % % September 8, 2003. Version 2.3. % http://www.cise.ufl.edu/research/sparse/colamd/ % %------------------------------------------------------------------------------- % Print the introduction, the help info, and compile the mexFunctions %------------------------------------------------------------------------------- more on fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Colamd/symamd demo.') ; fprintf (1, '\n-----------------------------------------------------------\n') ; help colamd_demo ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Colamd help information:') ; fprintf (1, '\n-----------------------------------------------------------\n') ; help colamd ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Symamd help information:') ; fprintf (1, '\n-----------------------------------------------------------\n') ; help symamd ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'colamd and symamd mexFunctions will now be compiled.') ; fprintf (1, '\n-----------------------------------------------------------\n') ; s = input ('\n\nHit any key (or type ''n'' to skip compilation): ', 's') ; if (~strcmp (s, 'n')) disp ('mex -O colamdmex.c colamd.c') ; mex -O colamdmex.c colamd.c disp ('mex -O symamdmex.c colamd.c') ; mex -O symamdmex.c colamd.c end %------------------------------------------------------------------------------- % Solving Ax=b %------------------------------------------------------------------------------- n = 100 ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Solving Ax=b for a small %d-by-%d random matrix:', n, n) ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, '\nNOTE: Random sparse matrices are AWFUL test cases.\n') ; fprintf (1, 'They''re just easy to generate in a demo.\n') ; % set up the system rand ('state', 0) ; randn ('state', 0) ; spparms ('default') ; A = sprandn (n, n, 5/n) + speye (n) ; b = (1:n)' ; fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colamd:\n') ; q = colamd (A) ; I = speye (n) ; Q = I (:, q) ; [L,U,P] = lu (A*Q) ; fl = luflops (L, U) ; x = Q * (U \ (L \ (P * b))) ; fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ; fprintf (1, 'residual: %e\n', norm (A*x-b)); fprintf (1, '\n\nSolving via lu (PAQ = LU), where Q is from colmmd:\n') ; q = colmmd (A) ; I = speye (n) ; Q = I (:, q) ; [L,U,P] = lu (A*Q) ; fl = luflops (L, U) ; x = Q * (U \ (L \ (P * b))) ; fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ; fprintf (1, 'residual: %e\n', norm (A*x-b)); fprintf (1, '\n\nSolving via lu (PA = LU), without regard for sparsity:\n') ; [L,U,P] = lu (A) ; fl = luflops (L, U) ; x = U \ (L \ (P * b)) ; fprintf (1, '\nFlop count for [L,U,P] = lu (A*Q): %d\n', fl) ; fprintf (1, 'residual: %e\n', norm (A*x-b)); %------------------------------------------------------------------------------- % Large demo for colamd %------------------------------------------------------------------------------- fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Large demo for colamd (symbolic analysis only):') ; fprintf (1, '\n-----------------------------------------------------------\n') ; rand ('state', 0) ; randn ('state', 0) ; spparms ('default') ; n = 1000 ; fprintf (1, 'Generating a random %d-by-%d sparse matrix.\n', n, n) ; A = sprandn (n, n, 5/n) + speye (n) ; fprintf (1, '\n\nUnordered matrix:\n') ; lnz = symbfact (A, 'col') ; fprintf (1, 'nz in Cholesky factors of A''A: %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A''A: %d\n', sum (lnz.^2)) ; tic ; p = colamd (A) ; t = toc ; lnz = symbfact (A (:,p), 'col') ; fprintf (1, '\n\nColamd run time: %f\n', t) ; fprintf (1, 'colamd ordering quality: \n') ; fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p): %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', sum (lnz.^2)) ; tic ; p = colmmd (A) ; t = toc ; lnz = symbfact (A (:,p), 'col') ; fprintf (1, '\n\nColmmd run time: %f\n', t) ; fprintf (1, 'colmmd ordering quality: \n') ; fprintf (1, 'nz in Cholesky factors of A(:,p)''A(:,p): %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A(:,p)''A(:,p): %d\n', sum (lnz.^2)) ; %------------------------------------------------------------------------------- % Large demo for symamd %------------------------------------------------------------------------------- fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Large demo for symamd (symbolic analysis only):') ; fprintf (1, '\n-----------------------------------------------------------\n') ; fprintf (1, 'Generating a random symmetric %d-by-%d sparse matrix.\n', n, n) ; A = A+A' ; fprintf (1, '\n\nUnordered matrix:\n') ; lnz = symbfact (A, 'sym') ; fprintf (1, 'nz in Cholesky factors of A: %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A: %d\n', sum (lnz.^2)) ; tic ; p = symamd (A) ; t = toc ; lnz = symbfact (A (p,p), 'sym') ; fprintf (1, '\n\nSymamd run time: %f\n', t) ; fprintf (1, 'symamd ordering quality: \n') ; fprintf (1, 'nz in Cholesky factors of A(p,p): %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ; tic ; p = symmmd (A) ; t = toc ; lnz = symbfact (A (p,p), 'sym') ; fprintf (1, '\n\nSymmmd run time: %f\n', t) ; fprintf (1, 'symmmd ordering quality: \n') ; fprintf (1, 'nz in Cholesky factors of A(p,p): %d\n', sum (lnz)) ; fprintf (1, 'flop count for Cholesky of A(p,p): %d\n', sum (lnz.^2)) ; more off SHAR_EOF fi # end of overwriting check if test -f 'colamd_test.m' then echo shar: will not over-write existing file "'colamd_test.m'" else cat << "SHAR_EOF" > 'colamd_test.m' function colamd_test (spumoni) % % colamd_test (spumoni) % % COLAMD and SYMAMD testing function. Here we try to give colamd and symamd % every possible type of matrix and erroneous input that they may encounter. % We want either a valid permutation returned or we want them to fail % gracefully. If the optional spumoni argument is present, additional % diagnostic messages are printed during the test: % % spumoni: % 0: very little, just a progress meter (default) % 1: lots % 2: far too much % % You are prompted as to whether or not the colamd and symand routines and % the test mexFunctions are to be compiled. % % Tim Davis % September 8, 2003. Version 2.3. % http://www.cise.ufl.edu/research/sparse/colamd/ help colamd_test s = input (... 'Compile colamd, symand, and the test codes? (y/n, default is yes): ', 's') ; do_compile = 1 ; if (~isempty (s)) if (s (1) == 'n' | s (1) == 'N') do_compile = 0 ; end end if (do_compile) fprintf ('Compiling colamd, symamd, and test mexFunctions.\n') ; fprintf ('mex -O colamdmex.c colamd.c\n') ; mex -O colamdmex.c colamd.c fprintf ('mex -O symamdmex.c colamd.c\n') ; mex -O symamdmex.c colamd.c fprintf ('mex -O colamdtestmex.c colamd.c\n') ; mex -O colamdtestmex.c colamd.c fprintf ('mex -O symamdtestmex.c colamd.c\n') ; mex -O symamdtestmex.c colamd.c fprintf ('Done compiling.\n') ; end fprintf ('\nThe following codes will be tested:\n') ; which colamd which symamd which colamdmex which symamdmex fprintf ('\nStarting the tests. Please be patient.\n') ; rand ('state', 0) ; randn ('state', 0) ; if (nargin < 1) spumoni = 0 ; end old = spparms ('spumoni') ; spparms ('spumoni', spumoni) ; fprintf ('Null matrices') ; A = zeros (0,0) ; A = sparse (A) ; [p, stats] = colamd (A, [.5 .5 spumoni]) ; check_perm (p, A) ; [p, stats] = symamd (A, [.5 spumoni]) ; check_perm (p, A) ; A = zeros (0, 100) ; A = sparse (A) ; [p, stats] = colamd (A, [.5 .5 spumoni]) ; check_perm (p, A) ; A = zeros (100, 0) ; A = sparse (A) ; [p, stats] = colamd (A, [.5 .5 spumoni]) ; check_perm (p, A) ; fprintf (' OK\n') ; fprintf ('Matrices with a few dense row/cols\n') ; for trial = 1:20 % random square unsymmetric matrix A = rand_matrix (1000, 1000, 1, 10, 20) ; [m n] = size (A) ; for tol = [0 0.001 0.005 0.01:0.01:0.10 0.5:.1:1] [p, stats] = colamd (A, [tol tol spumoni]) ; check_perm (p, A) ; B = A + A' ; [p, stats] = symamd (B, [tol spumoni]) ; check_perm (p, A) ; [p, stats] = colamd (A, [tol 1 spumoni]) ; check_perm (p, A) ; [p, stats] = colamd (A, [1 tol spumoni]) ; check_perm (p, A) ; fprintf ('.') ; end end fprintf (' OK\n') ; fprintf ('General matrices\n') ; for trial = 1:400 % matrix of random mtype mtype = irand (3) ; A = rand_matrix (2000, 2000, mtype, 0, 0) ; [m n] = size (A) ; p = colamd (A) ; check_perm (p, A) ; if (mtype == 3) p = symamd (A) ; check_perm (p, A) ; end fprintf ('.') ; end fprintf (' OK\n') ; fprintf ('Test error handling with invalid inputs\n') ; % Check different erroneous input. for trial = 1:30 A = rand_matrix (1000, 1000, 2, 0, 0) ; [m n] = size (A) ; for err = 1:13 p = Tcolamd (A, [1 1 spumoni 0 err]) ; if p ~= -1 check_perm (p, A) ; end if (err == 1) % check different (valid) input args to colamd p = Acolamd (A, spumoni) ; p2 = Acolamd (A, spumoni, [.5 .5 spumoni 0 0]) ; if (any (p ~= p2)) error ('colamd: mismatch 1!') ; end [p2 stats] = Acolamd (A, spumoni) ; if (any (p ~= p2)) error ('colamd: mismatch 2!') ; end [p2 stats] = Acolamd (A, spumoni, [.5 .5 spumoni 0 0]) ; if (any (p ~= p2)) error ('colamd: mismatch 3!') ; end end B = A'*A ; p = Tsymamd (B, [1 spumoni err]) ; if p ~= -1 check_perm (p, A) ; end if (err == 1) % check different (valid) input args to symamd p = Asymamd (B, spumoni) ; check_perm (p, A) ; p2 = Asymamd (B, spumoni, [.5 spumoni 0]) ; if (any (p ~= p2)) error ('symamd: mismatch 1!') ; end [p2 stats] = Asymamd (B, spumoni) ; if (any (p ~= p2)) error ('symamd: mismatch 2!') ; end [p2 stats] = Asymamd (B, spumoni, [.5 spumoni 0]) ; if (any (p ~= p2)) error ('symamd: mismatch 3!') ; end end fprintf ('.') ; end end fprintf (' OK\n') ; fprintf ('Matrices with a few empty columns\n') ; for trial = 1:400 % some are square, some are rectangular n = 0 ; while (n < 5) A = rand_matrix (1000, 1000, irand (2), 0, 0) ; [m n] = size (A) ; end % Add 5 null columns at random locations. null_col = randperm (n) ; null_col = sort (null_col (1:5)) ; A (:, null_col) = 0 ; % Order the matrix and make sure that the null columns are ordered last. [p, stats] = colamd (A, [1 1 spumoni]) ; check_perm (p, A) ; if (stats (2) ~= 5) error ('colamd: wrong number of null columns') ; end if (any (null_col ~= p ((n-4):n))) error ('colamd: Null cols are not ordered last in natural order') ; end fprintf ('.') ; end fprintf (' OK\n') ; fprintf ('Matrices with a few empty rows and columns\n') ; for trial = 1:400 % symmetric matrices n = 0 ; while (n < 5) A = rand_matrix (1000, 1000, 3, 0, 0) ; [m n] = size (A) ; end % Add 5 null columns and rows at random locations. null_col = randperm (n) ; null_col = sort (null_col (1:5)) ; A (:, null_col) = 0 ; A (null_col, :) = 0 ; % Order the matrix and make sure that the null rows/cols are ordered last. [p,stats] = symamd (A, [1 spumoni]) ; check_perm (p, A) ; % find actual number of null rows and columns Alo = tril (A, -1) ; nnull = length (find (sum (Alo') == 0 & sum (Alo) == 0)) ; if (stats (2) ~= nnull | nnull < 5) error ('symamd: wrong number of null columns') ; end if (any (null_col ~= p ((n-4):n))) error ('symamd: Null cols are not ordered last in natural order') ; end fprintf ('.') ; end fprintf (' OK\n') ; fprintf ('Matrices with a few empty rows\n') ; % Test matrices with null rows inserted. for trial = 1:400 m = 0 ; while (m < 5) A = rand_matrix (1000, 1000, 2, 0, 0) ; [m n] = size (A) ; end % Add 5 null rows at random locations. null_row = randperm (m) ; null_row = sort (null_row (1:5)) ; A (null_row, :) = 0 ; p = colamd (A, [.5 .5 spumoni]) ; check_perm (p, A) ; if (stats (1) ~= 5) error ('colamd: wrong number of null rows') ; end fprintf ('.') ; end fprintf (' OK\n') ; fprintf ('\ncolamd and symamd: all tests passed\n\n') ; spparms ('spumoni', old) ; %------------------------------------------------------------------------------- % Acolamd: compare colamd and Tcolamd results %------------------------------------------------------------------------------- function [p,stats] = Acolamd (S, spumoni, knobs) if (nargin < 2) spumoni = 0 ; end if (nargin < 3) if (nargout == 1) [p] = colamd (S) ; [p1] = Tcolamd (S, [.5 .5 spumoni 0 0]) ; else [p, stats] = colamd (S) ; [p1, stats1] = Tcolamd (S, [.5 .5 spumoni 0 0]) ; end else if (nargout == 1) [p] = colamd (S, knobs (1:3)) ; [p1] = Tcolamd (S, knobs) ; else [p, stats] = colamd (S, knobs (1:3)) ; [p1, stats1] = Tcolamd (S, knobs) ; end end check_perm (p, S) ; check_perm (p1, S) ; if (any (p1 ~= p)) error ('Acolamd mismatch!') ; end %------------------------------------------------------------------------------- % Asymamd: compare symamd and Tsymamd results %------------------------------------------------------------------------------- function [p,stats] = Asymamd (S, spumoni, knobs) if (nargin < 2) spumoni = 0 ; end if (nargin < 3) if (nargout == 1) [p] = symamd (S) ; [p1] = Tsymamd (S, [.5 spumoni 0]) ; else [p, stats] = symamd (S) ; [p1, stats1] = Tsymamd (S, [.5 spumoni 0]) ; end else if (nargout == 1) [p] = symamd (S, knobs (1:2)) ; [p1] = Tsymamd (S, knobs) ; else [p, stats] = symamd (S, knobs (1:2)) ; [p1, stats1] = Tsymamd (S, knobs) ; end end if (any (p1 ~= p)) error ('Asymamd mismatch!') ; end %------------------------------------------------------------------------------- % check_perm: check for a valid permutation vector %------------------------------------------------------------------------------- function check_perm (p, A) if (isempty (A) & isempty (p)) % empty permutation vectors of empty matrices are OK return end if (isempty (p)) error ('bad permutation: cannot be empty') ; end [m n] = size (A) ; [pm pn] = size (p) ; if (pn == 1) % force p to be a row vector p = p' ; [pm pn] = size (p) ; end if (n ~= pn) error ('bad permutation: wrong size') ; end if (pm ~= 1) ; % p must be a vector error ('bad permutation: not a vector') ; else if (any (sort (p) - (1:pn))) error ('bad permutation') ; end end %------------------------------------------------------------------------------- % irand: return a random integer between 1 and n %------------------------------------------------------------------------------- function i = irand (n) i = min (n, 1 + floor (rand * n)) ; %------------------------------------------------------------------------------- % rand_matrix: return a random sparse matrix %------------------------------------------------------------------------------- function A = rand_matrix (nmax, mmax, mtype, drows, dcols) % % A = rand_matrix (nmax, mmax, mtype, drows, dcols) % % A binary matrix of random size, at most nmax-by-mmax, with drows dense rows % and dcols dense columns. % % mtype 1: square unsymmetric (mmax is ignored) % mtype 2: rectangular % mtype 3: symmetric (mmax is ignored) n = irand (nmax) ; if (mtype ~= 2) % square m = n ; else m = irand (mmax) ; end A = sprand (m, n, 10 / max (m,n)) ; if (drows > 0) % add dense rows for k = 1:drows i = irand (m) ; nz = irand (n) ; p = randperm (n) ; p = p (1:nz) ; A (i,p) = 1 ; end end if (dcols > 0) % add dense cols for k = 1:dcols j = irand (n) ; nz = irand (m) ; p = randperm (m) ; p = p (1:nz) ; A (p,j) = 1 ; end end A = spones (A) ; % ensure that there are no empty columns d = find (full (sum (A)) == 0) ; A (m,d) = 1 ; % ensure that there are no empty rows d = find (full (sum (A,2)) == 0) ; A (d,n) = 1 ; if (mtype == 3) % symmetric A = A + A' + speye (n) ; end A = spones (A) ; %------------------------------------------------------------------------------- % Tcolamd: run colamd in a testing mode %------------------------------------------------------------------------------- function [p,stats] = Tcolamd (S, knobs) if (nargout <= 1 & nargin == 1) p = colamdtestmex (S) ; elseif (nargout <= 1 & nargin == 2) p = colamdtestmex (S, knobs) ; elseif (nargout == 2 & nargin == 1) [p, stats] = colamdtestmex (S) ; elseif (nargout == 2 & nargin == 2) [p, stats] = colamdtestmex (S, knobs) ; else error ('colamd: incorrect number of input and/or output arguments') ; end if (p (1) ~= -1) [ignore, q] = sparsfun ('coletree', S (:,p)) ; p = p (q) ; check_perm (p, S) ; end %------------------------------------------------------------------------------- % Tsymamd: run symamd in a testing mode %------------------------------------------------------------------------------- function [p, stats] = Tsymamd (S, knobs) if (nargout <= 1 & nargin == 1) p = symamdtestmex (S) ; elseif (nargout <= 1 & nargin == 2) p = symamdtestmex (S, knobs) ; elseif (nargout == 2 & nargin == 1) [p, stats] = symamdtestmex (S) ; elseif (nargout == 2 & nargin == 2) [p, stats] = symamdtestmex (S, knobs) ; else error ('symamd: incorrect number of input and/or output arguments') ; end if (p (1) ~= -1) [ignore, q] = sparsfun ('symetree', S (p,p)) ; p = p (q) ; check_perm (p, S) ; end SHAR_EOF fi # end of overwriting check if test -f 'colamdtestmex.c' then echo shar: will not over-write existing file "'colamdtestmex.c'" else cat << "SHAR_EOF" > 'colamdtestmex.c' /* ========================================================================== */ /* === colamdtest mexFunction =============================================== */ /* ========================================================================== */ /* This MATLAB mexFunction is for testing only. It is not meant for production use. See colamdmex.c instead. Usage: [ P, stats ] = colamdtest (A, knobs) ; Returns a permutation vector P such that the LU factorization of A (:,P) tends to be sparser than that of A. The Cholesky factorization of (A (:,P))'*(A (:,P)) will also tend to be sparser than that of A'*A. This routine provides the same functionality as COLMMD, but is much faster and returns a better permutation vector. Note that the COLMMD m-file in MATLAB 5.2 also performs a column elimination tree post-ordering. This mexFunction does not do this post-ordering. This mexFunction is a replacement for the p = sparsfun ('colmmd', A) operation. The knobs and stats vectors are optional: knobs (1) rows with more than (knobs (1))*n_col entries are removed prior to ordering. If knobs is not present, then the default is used (0.5). knobs (2) columns with more than (knobs (2))*n_row entries are removed prior to ordering, and placed last in the column permutation. If knobs is not present, then the default is used (0.5). knobs (3) print level, similar to spparms ('spumoni') knobs (4) for testing only. Controls the workspace used by colamd. knobs (5) for testing only. Controls how the input matrix is jumbled prior to calling colamd, to test its error handling capability. stats (1) the number of dense (or empty) rows ignored stats (2) the number of dense (or empty) columms. These are ordered last, in their natural order. stats (3) the number of garbage collections performed. stats (4) return status: 0: matrix is a valid MATLAB matrix. 1: matrix has duplicate entries or unsorted columns. This should not occur in a valid MATLAB matrix, but the ordering proceeded anyway by sorting the row indices in each column and by ignoring the duplicate row indices in each column. See stats (5:7) for more information. stats (5) highest numbered column that is unsorted or has duplicate entries (zero if none) stats (6) last seen duplicate or unsorted row index (zero if none) stats (7) number of duplicate or unsorted row indices Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Notice: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c file) for the License. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/colamdtestmex.c file. It requires the colamd.c and colamd.h files. */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "colamd.h" #include "mex.h" #include "matrix.h" #include #include static void dump_matrix ( int A [ ], int p [ ], int n_row, int n_col, int Alen, int limit ) ; /* ========================================================================== */ /* === colamd mexFunction =================================================== */ /* ========================================================================== */ void mexFunction ( /* === Parameters ======================================================= */ int nlhs, /* number of left-hand sides */ mxArray *plhs [], /* left-hand side matrices */ int nrhs, /* number of right--hand sides */ const mxArray *prhs [] /* right-hand side matrices */ ) { /* === Local variables ================================================== */ int *A ; /* colamd's copy of the matrix, and workspace */ int *p ; /* colamd's copy of the column pointers */ int Alen ; /* size of A */ int n_col ; /* number of columns of A */ int n_row ; /* number of rows of A */ int nnz ; /* number of entries in A */ int full ; /* TRUE if input matrix full, FALSE if sparse */ double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */ double *out_perm ; /* output permutation vector */ double *out_stats ; /* output stats vector */ double *in_knobs ; /* input knobs vector */ int i ; /* loop counter */ mxArray *Ainput ; /* input matrix handle */ int spumoni ; /* verbosity variable */ int stats2 [COLAMD_STATS] ; /* stats for colamd */ int *cp, *cp_end, result, col, length ; int *stats ; stats = stats2 ; /* === Check inputs ===================================================== */ if (nrhs < 1 || nrhs > 2 || nlhs < 0 || nlhs > 2) { mexErrMsgTxt ( "colamd: incorrect number of input and/or output arguments") ; } if (nrhs != 2) { mexErrMsgTxt ("colamdtest: knobs are required") ; } /* for testing we require all 5 knobs */ if (mxGetNumberOfElements (prhs [1]) < 5) { mexErrMsgTxt ("colamd: must have all 5 knobs for testing") ; } /* === Get knobs ======================================================== */ colamd_set_defaults (knobs) ; spumoni = 0 ; /* check for user-passed knobs */ if (nrhs == 2) { in_knobs = mxGetPr (prhs [1]) ; i = mxGetNumberOfElements (prhs [1]) ; if (i > 0) knobs [COLAMD_DENSE_ROW] = in_knobs [COLAMD_DENSE_ROW] ; if (i > 1) knobs [COLAMD_DENSE_COL] = in_knobs [COLAMD_DENSE_COL] ; if (i > 2) spumoni = (int) in_knobs [2] ; } /* print knob settings if spumoni is set */ if (spumoni > 0) { mexPrintf ("colamd: dense row fraction: %f\n", knobs [COLAMD_DENSE_ROW]) ; mexPrintf ("colamd: dense col fraction: %f\n", knobs [COLAMD_DENSE_COL]) ; } /* === If A is full, convert to a sparse matrix ========================= */ Ainput = (mxArray *) prhs [0] ; if (mxGetNumberOfDimensions (Ainput) != 2) { mexErrMsgTxt ("colamd: input matrix must be 2-dimensional") ; } full = !mxIsSparse (Ainput) ; if (full) { mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ; } /* === Allocate workspace for colamd ==================================== */ /* get size of matrix */ n_row = mxGetM (Ainput) ; n_col = mxGetN (Ainput) ; /* get column pointer vector so we can find nnz */ p = (int *) mxCalloc (n_col+1, sizeof (int)) ; (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (int)) ; nnz = p [n_col] ; Alen = colamd_recommended (nnz, n_row, n_col) ; /* === Modify size of Alen if testing ======================================= */ /* knobs [3] amount of workspace given to colamd. < 0 : TIGHT memory > 0 : MIN + knob [3] - 1 == 0 : RECOMMENDED memory */ /* Here only for testing */ #ifdef MIN #undef MIN #endif #define MIN(a,b) (((a) < (b)) ? (a) : (b)) #define COLAMD_MIN_MEMORY(nnz,n_row,n_col) \ (2 * (nnz) + COLAMD_C (n_col) + COLAMD_R (n_row)) /* get knob [3], if negative */ if (in_knobs [3] < 0) { Alen = COLAMD_MIN_MEMORY (nnz, n_row, n_col) + n_col ; } else if (in_knobs [3] > 0) { Alen = COLAMD_MIN_MEMORY (nnz, n_row, n_col) + in_knobs [3] - 1 ; } /* otherwise, we use the recommended amount set above */ /* === Copy input matrix into workspace ================================= */ A = (int *) mxCalloc (Alen, sizeof (int)) ; (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (int)) ; if (full) { mxDestroyArray (Ainput) ; } /* === Jumble matrix ======================================================== */ /* knobs [4] FOR TESTING ONLY: Specifies how to jumble matrix 0 : No jumbling 1 : Make n_row less than zero 2 : Make first pointer non-zero 3 : Make column pointers not non-decreasing 4 : Make a column pointer greater or equal to Alen 5 : Make row indices not strictly increasing 6 : Make a row index greater or equal to n_row 7 : Set A = NULL 8 : Set p = NULL 9 : Repeat row index 10: make row indices not sorted 11: jumble columns massively (note this changes the pattern of the matrix A.) 12: Set stats = NULL 13: Make n_col less than zero */ /* jumble appropriately */ switch ((int) in_knobs [4]) { case 0 : if (spumoni > 0) { mexPrintf ("colamdtest: no errors expected\n") ; } result = 1 ; /* no errors */ break ; case 1 : if (spumoni > 0) { mexPrintf ("colamdtest: nrow out of range\n") ; } result = 0 ; /* nrow out of range */ n_row = -1 ; break ; case 2 : if (spumoni > 0) { mexPrintf ("colamdtest: p [0] nonzero\n") ; } result = 0 ; /* p [0] must be zero */ p [0] = 1 ; break ; case 3 : if (spumoni > 0) { mexPrintf ("colamdtest: negative length last column\n") ; } result = (n_col == 0) ; /* p must be monotonically inc. */ p [n_col] = p [0] ; break ; case 4 : if (spumoni > 0) { mexPrintf ("colamdtest: Alen too small\n") ; } result = 0 ; /* out of memory */ p [n_col] = Alen ; break ; case 5 : if (spumoni > 0) { mexPrintf ("colamdtest: row index out of range (-1)\n") ; } if (nnz > 0) /* row index out of range */ { result = 0 ; A [nnz-1] = -1 ; } else { if (spumoni > 0) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 6 : if (spumoni > 0) { mexPrintf ("colamdtest: row index out of range (n_row)\n") ; } if (nnz > 0) /* row index out of range */ { if (spumoni > 0) { mexPrintf ("Changing A[nnz-1] from %d to %d\n", A [nnz-1], n_row) ; } result = 0 ; A [nnz-1] = n_row ; } else { if (spumoni > 0) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 7 : if (spumoni > 0) { mexPrintf ("colamdtest: A not present\n") ; } result = 0 ; /* A not present */ A = (int *) NULL ; break ; case 8 : if (spumoni > 0) { mexPrintf ("colamdtest: p not present\n") ; } result = 0 ; /* p not present */ p = (int *) NULL ; break ; case 9 : if (spumoni > 0) { mexPrintf ("colamdtest: duplicate row index\n") ; } result = 1 ; /* duplicate row index */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { A [p [col]] = A [p [col] + 1] ; if (spumoni > 0) { mexPrintf ("Made duplicate row %d in col %d\n", A [p [col] + 1], col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, col+2) ; } break ; case 10 : if (spumoni > 0) { mexPrintf ("colamdtest: unsorted column\n") ; } result = 1 ; /* jumbled columns */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { i = A[p [col]] ; A [p [col]] = A[p [col] + 1] ; A [p [col] + 1] = i ; if (spumoni > 0) { mexPrintf ("Unsorted column %d \n", col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, col+2) ; } break ; case 11 : if (spumoni > 0) { mexPrintf ("colamdtest: massive jumbling\n") ; } result = 1 ; /* massive jumbling, but no errors */ srand (1) ; for (i = 0 ; i < n_col ; i++) { cp = &A [p [i]] ; cp_end = &A [p [i+1]] ; while (cp < cp_end) { *cp++ = rand() % n_row ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, Alen, n_col) ; } break ; case 12 : if (spumoni > 0) { mexPrintf ("colamdtest: stats not present\n") ; } result = 0 ; /* stats not present */ stats = (int *) NULL ; break ; case 13 : if (spumoni > 0) { mexPrintf ("colamdtest: ncol out of range\n") ; } result = 0 ; /* ncol out of range */ n_col = -1 ; break ; } /* === Order the columns (destroys A) =================================== */ if (!colamd (n_row, n_col, Alen, A, p, knobs, stats)) { /* return p = -1 if colamd failed */ plhs [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; out_perm [0] = -1 ; mxFree (p) ; mxFree (A) ; if (spumoni > 0 || result) { colamd_report (stats) ; } if (result) { mexErrMsgTxt ("colamd should have returned TRUE\n") ; } return ; /* mexErrMsgTxt ("colamd error!") ; */ } if (!result) { colamd_report (stats) ; mexErrMsgTxt ("colamd should have returned FALSE\n") ; } mxFree (A) ; /* === Return the permutation vector ==================================== */ plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; for (i = 0 ; i < n_col ; i++) { /* colamd is 0-based, but MATLAB expects this to be 1-based */ out_perm [i] = p [i] + 1 ; } mxFree (p) ; /* === Return the stats vector ========================================== */ /* print stats if spumoni > 0 */ if (spumoni > 0) { colamd_report (stats) ; } if (nlhs == 2) { plhs [1] = mxCreateDoubleMatrix (1, COLAMD_STATS, mxREAL) ; out_stats = mxGetPr (plhs [1]) ; for (i = 0 ; i < COLAMD_STATS ; i++) { out_stats [i] = stats [i] ; } /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ /* note that this correction doesn't occur if symamd returns FALSE */ out_stats [COLAMD_INFO1] ++ ; out_stats [COLAMD_INFO2] ++ ; } } static void dump_matrix ( int A [ ], int p [ ], int n_row, int n_col, int Alen, int limit ) { int col, k, row ; mexPrintf ("dump matrix: nrow %d ncol %d Alen %d\n", n_row, n_col, Alen) ; for (col = 0 ; col < MIN (n_col, limit) ; col++) { mexPrintf ("column %d, p[col] %d, p [col+1] %d, length %d\n", col, p [col], p [col+1], p [col+1] - p [col]) ; for (k = p [col] ; k < p [col+1] ; k++) { row = A [k] ; mexPrintf (" %d", row) ; } mexPrintf ("\n") ; } } SHAR_EOF fi # end of overwriting check if test -f 'luflops.m' then echo shar: will not over-write existing file "'luflops.m'" else cat << "SHAR_EOF" > 'luflops.m' function fl = luflops (L, U) % % fl = luflops (L,U) % % Given a sparse LU factorization (L and U), return the flop count required % by a conventional LU factorization algorithm to compute it. L and U can % be either sparse or full matrices. L must be lower triangular and U must % be upper triangular. Do not attempt to use this on the permuted L from % [L,U] = lu (A). Instead, use [L,U,P] = lu (A) or [L,U,P,Q] = lu (A). % % Note that there is a subtle undercount in this estimate. Suppose A is % completely dense, but during LU factorization exact cancellation occurs, % causing some of the entries in L and U to become identically zero. The % flop count returned by this routine is an undercount. There is a simple % way to fix this (L = spones (L) + spones (tril (A))), but the fix is partial. % It can also occur that some entry in L is a "symbolic" fill-in (zero in % A, but a fill-in entry and thus must be computed), but numerically % zero. The only way to get a reliable LU factorization would be to do a % purely symbolic factorization of A. This cannot be done with % symbfact (A, 'col'). % % See NA Digest, Vol 00, #50, Tuesday, Dec. 5, 2000 % % Tim Davis, Sept. 23, 2002. Written for MATLAB 6.5. Lnz = full (sum (spones (L))) - 1 ; % off diagonal nz in cols of L Unz = full (sum (spones (U')))' - 1 ; % off diagonal nz in rows of U fl = 2*Lnz*Unz + sum (Lnz) ; SHAR_EOF fi # end of overwriting check if test -f 'symamdtestmex.c' then echo shar: will not over-write existing file "'symamdtestmex.c'" else cat << "SHAR_EOF" > 'symamdtestmex.c' /* ========================================================================== */ /* === symamdtest mexFunction =============================================== */ /* ========================================================================== */ /* This MATLAB mexFunction is for testing only. It is not meant for production use. See symamdmex.c instead. Usage: [ P, stats ] = symamdtest (A, knobs) ; Returns a permutation vector P such that the Cholesky factorization of A (P,P) tends to be sparser than that of A. This routine provides the same functionality as SYMMMD, but tends to be much faster and tends to return a better permutation vector. Note that the SYMMMD m-file in MATLAB 5.2 also performs a symmetric elimination tree post-ordering. This mexFunction does not do this post-ordering. This mexFunction is a replacement for the p = sparsfun ('symmmd', A) operation. A must be square, and is assummed to have a symmetric nonzero pattern. Only the nonzero pattern of the lower triangular portion of A is accessed. This routine constructs a matrix M such that the nonzero pattern of M'M is equal to A (assuming A has symmetric pattern), and then performs a column ordering of M using colamd. The knobs and stats vectors are optional: knobs (1) rows and columns with more than (knobs(1))*n entries are removed prior to ordering, and placed last in the output ordering. If knobs is not present, then the default of 0.5 is used. knobs (2) print level, similar to spparms ('spumoni') knobs (3) for testing only. Controls how the input matrix is jumbled prior to calling symamd, to test its error handling capability. stats (1) the number of dense (or empty) rows and columms. These are ordered last, in their natural order. stats (2) (same as stats (1)) stats (3) the number of garbage collections performed. stats (4) return status: 0: matrix is a valid MATLAB matrix. 1: matrix has duplicate entries or unsorted columns. This should not occur in a valid MATLAB matrix, but the ordering proceeded anyway by ignoring the duplicate row indices in each column. See stats (5:7) for more information. stats (5) highest numbered column that is unsorted or has duplicate entries (zero if none) stats (6) last seen duplicate or unsorted row index (zero if none) stats (7) number of duplicate or unsorted row indices Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Notice: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c file) for the License. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/symamdtestmex.c file. It requires the colamd.c and colamd.h files. */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "colamd.h" #include "mex.h" #include "matrix.h" #include #include static void dump_matrix ( int A [ ], int p [ ], int n_row, int n_col, int Alen, int limit ) ; /* ========================================================================== */ /* === symamd mexFunction =================================================== */ /* ========================================================================== */ void mexFunction ( /* === Parameters ======================================================= */ int nlhs, /* number of left-hand sides */ mxArray *plhs [], /* left-hand side matrices */ int nrhs, /* number of right--hand sides */ const mxArray *prhs [] /* right-hand side matrices */ ) { /* === Local variables ================================================== */ int *perm ; /* column ordering of M and ordering of A */ int *A ; /* row indices of input matrix A */ int *p ; /* column pointers of input matrix A */ int n_col ; /* number of columns of A */ int n_row ; /* number of rows of A */ int full ; /* TRUE if input matrix full, FALSE if sparse */ double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */ double *out_perm ; /* output permutation vector */ double *out_stats ; /* output stats vector */ double *in_knobs ; /* input knobs vector */ int i ; /* loop counter */ mxArray *Ainput ; /* input matrix handle */ int spumoni ; /* verbosity variable */ int stats2 [COLAMD_STATS] ; /* stats for symamd */ int *cp, *cp_end, result, nnz, col, length ; int *stats ; stats = stats2 ; /* === Check inputs ===================================================== */ if (nrhs < 1 || nrhs > 2 || nlhs < 0 || nlhs > 2) { mexErrMsgTxt ( "symamd: incorrect number of input and/or output arguments.") ; } if (nrhs != 2) { mexErrMsgTxt ("symamdtest: knobs are required") ; } /* for testing we require all 3 knobs */ if (mxGetNumberOfElements (prhs [1]) < 3) { mexErrMsgTxt ("symamdtest: must have all 3 knobs for testing") ; } /* === Get knobs ======================================================== */ colamd_set_defaults (knobs) ; spumoni = 0 ; /* check for user-passed knobs */ if (nrhs == 2) { in_knobs = mxGetPr (prhs [1]) ; i = mxGetNumberOfElements (prhs [1]) ; if (i > 0) knobs [COLAMD_DENSE_ROW] = in_knobs [COLAMD_DENSE_ROW] ; if (i > 1) spumoni = (int) in_knobs [1] ; } /* print knob settings if spumoni > 0 */ if (spumoni > 0) { mexPrintf ("symamd: dense row/col fraction: %f\n", knobs [COLAMD_DENSE_ROW]) ; } /* === If A is full, convert to a sparse matrix ========================= */ Ainput = (mxArray *) prhs [0] ; if (mxGetNumberOfDimensions (Ainput) != 2) { mexErrMsgTxt ("symamd: input matrix must be 2-dimensional.") ; } full = !mxIsSparse (Ainput) ; if (full) { mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ; } /* === Allocate workspace for symamd ==================================== */ /* get size of matrix */ n_row = mxGetM (Ainput) ; n_col = mxGetN (Ainput) ; if (n_col != n_row) { mexErrMsgTxt ("symamd: matrix must be square.") ; } /* p = mxGetJc (Ainput) ; */ p = (int *) mxCalloc (n_col+1, sizeof (int)) ; (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (int)) ; nnz = p [n_col] ; if (spumoni > 0) { mexPrintf ("symamdtest: nnz %d\n", nnz) ; } /* A = mxGetIr (Ainput) ; */ A = (int *) mxCalloc (nnz+1, sizeof (int)) ; (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (int)) ; perm = (int *) mxCalloc (n_col+1, sizeof (int)) ; /* === Jumble matrix ======================================================== */ /* knobs [2] FOR TESTING ONLY: Specifies how to jumble matrix 0 : No jumbling 1 : (no errors) 2 : Make first pointer non-zero 3 : Make column pointers not non-decreasing 4 : (no errors) 5 : Make row indices not strictly increasing 6 : Make a row index greater or equal to n_row 7 : Set A = NULL 8 : Set p = NULL 9 : Repeat row index 10: make row indices not sorted 11: jumble columns massively (note this changes the pattern of the matrix A.) 12: Set stats = NULL 13: Make n_col less than zero */ /* jumble appropriately */ switch ((int) in_knobs [2]) { case 0 : if (spumoni > 0) { mexPrintf ("symamdtest: no errors expected\n") ; } result = 1 ; /* no errors */ break ; case 1 : if (spumoni > 0) { mexPrintf ("symamdtest: no errors expected (1)\n") ; } result = 1 ; break ; case 2 : if (spumoni > 0) { mexPrintf ("symamdtest: p [0] nonzero\n") ; } result = 0 ; /* p [0] must be zero */ p [0] = 1 ; break ; case 3 : if (spumoni > 0) { mexPrintf ("symamdtest: negative length last column\n") ; } result = (n_col == 0) ; /* p must be monotonically inc. */ p [n_col] = p [0] ; break ; case 4 : if (spumoni > 0) { mexPrintf ("symamdtest: no errors expected (4)\n") ; } result = 1 ; break ; case 5 : if (spumoni > 0) { mexPrintf ("symamdtest: row index out of range (-1)\n") ; } if (nnz > 0) /* row index out of range */ { result = 0 ; A [nnz-1] = -1 ; } else { if (spumoni > 0) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 6 : if (spumoni > 0) { mexPrintf ("symamdtest: row index out of range (ncol)\n") ; } if (nnz > 0) /* row index out of range */ { result = 0 ; A [nnz-1] = n_col ; } else { if (spumoni > 0) { mexPrintf ("Note: no row indices to put out of range\n") ; } result = 1 ; } break ; case 7 : if (spumoni > 0) { mexPrintf ("symamdtest: A not present\n") ; } result = 0 ; /* A not present */ A = (int *) NULL ; break ; case 8 : if (spumoni > 0) { mexPrintf ("symamdtest: p not present\n") ; } result = 0 ; /* p not present */ p = (int *) NULL ; break ; case 9 : if (spumoni > 0) { mexPrintf ("symamdtest: duplicate row index\n") ; } result = 1 ; /* duplicate row index */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { A [p [col+1]-2] = A [p [col+1] - 1] ; if (spumoni > 0) { mexPrintf ("Made duplicate row %d in col %d\n", A [p [col+1] - 1], col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, nnz, col+2) ; } break ; case 10 : if (spumoni > 0) { mexPrintf ("symamdtest: unsorted column\n") ; } result = 1 ; /* jumbled columns */ for (col = 0 ; col < n_col ; col++) { length = p [col+1] - p [col] ; if (length > 1) { i = A[p [col]] ; A [p [col]] = A[p [col] + 1] ; A [p [col] + 1] = i ; if (spumoni > 0) { mexPrintf ("Unsorted column %d \n", col) ; } break ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, nnz, col+2) ; } break ; case 11 : if (spumoni > 0) { mexPrintf ("symamdtest: massive jumbling\n") ; } result = 1 ; /* massive jumbling, but no errors */ srand (1) ; for (i = 0 ; i < n_col ; i++) { cp = &A [p [i]] ; cp_end = &A [p [i+1]] ; while (cp < cp_end) { *cp++ = rand() % n_row ; } } if (spumoni > 1) { dump_matrix (A, p, n_row, n_col, nnz, n_col) ; } break ; case 12 : if (spumoni > 0) { mexPrintf ("symamdtest: stats not present\n") ; } result = 0 ; /* stats not present */ stats = (int *) NULL ; break ; case 13 : if (spumoni > 0) { mexPrintf ("symamdtest: ncol out of range\n") ; } result = 0 ; /* ncol out of range */ n_col = -1 ; break ; } /* === Order the rows and columns of A (does not destroy A) ============= */ if (!symamd (n_col, A, p, perm, knobs, stats, &mxCalloc, &mxFree)) { /* return p = -1 if colamd failed */ plhs [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; out_perm [0] = -1 ; mxFree (p) ; mxFree (A) ; if (spumoni > 0 || result) { symamd_report (stats) ; } if (result) { mexErrMsgTxt ("symamd should have returned TRUE\n") ; } return ; /* mexErrMsgTxt ("symamd error!") ; */ } if (!result) { symamd_report (stats) ; mexErrMsgTxt ("symamd should have returned FALSE\n") ; } if (full) { mxDestroyArray (Ainput) ; } /* === Return the permutation vector ==================================== */ plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; for (i = 0 ; i < n_col ; i++) { /* symamd is 0-based, but MATLAB expects this to be 1-based */ out_perm [i] = perm [i] + 1 ; } mxFree (perm) ; /* === Return the stats vector ========================================== */ /* print stats if spumoni > 0 */ if (spumoni > 0) { symamd_report (stats) ; } if (nlhs == 2) { plhs [1] = mxCreateDoubleMatrix (1, COLAMD_STATS, mxREAL) ; out_stats = mxGetPr (plhs [1]) ; for (i = 0 ; i < COLAMD_STATS ; i++) { out_stats [i] = stats [i] ; } /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ /* note that this correction doesn't occur if symamd returns FALSE */ out_stats [COLAMD_INFO1] ++ ; out_stats [COLAMD_INFO2] ++ ; } } #ifdef MIN #undef MIN #endif #define MIN(a,b) (((a) < (b)) ? (a) : (b)) static void dump_matrix ( int A [ ], int p [ ], int n_row, int n_col, int Alen, int limit ) { int col, k, row ; mexPrintf ("dump matrix: nrow %d ncol %d Alen %d\n", n_row, n_col, Alen) ; if (!A) { mexPrintf ("A not present\n") ; return ; } if (!p) { mexPrintf ("p not present\n") ; return ; } for (col = 0 ; col < MIN (n_col, limit) ; col++) { mexPrintf ("column %d, p[col] %d, p [col+1] %d, length %d\n", col, p [col], p [col+1], p [col+1] - p [col]) ; for (k = p [col] ; k < p [col+1] ; k++) { row = A [k] ; mexPrintf (" %d", row) ; } mexPrintf ("\n") ; } } SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'colamd.m' then echo shar: will not over-write existing file "'colamd.m'" else cat << "SHAR_EOF" > 'colamd.m' function [p,stats] = colamd (S, knobs) %COLAMD Column approximate minimum degree permutation. % P = COLAMD (S) returns the column approximate minimum degree permutation % vector for the sparse matrix S. For a non-symmetric matrix S, S (:,P) % tends to have sparser LU factors than S. The Cholesky factorization of % S (:,P)' * S (:,P) also tends to be sparser than that of S'*S. COLAMD % tends to be faster than COLMMD and tends to return a better ordering. % % See also COLMMD, COLPERM, SPPARMS, SYMAMD, SYMMMD, SYMRCM. % % Usage: P = colamd (S) % P = colamd (S, knobs) % [P, stats] = colamd (S) % [P, stats] = colamd (S, knobs) % % knobs is an optional two-element input vector. If S is m-by-n, then % rows with more than (knobs (1))*n entries are ignored. Columns with % more than (knobs (2))*m entries are removed prior to ordering, and % ordered last in the output permutation P. If the knobs parameter is not % present, then 0.5 is used instead, for both knobs (1) and knobs (2). % knobs (3) controls the printing of statistics and error messages. % % stats is an optional 20-element output vector that provides data about the % ordering and the validity of the input matrix S. Ordering statistics are % in stats (1:3). stats (1) and stats (2) are the number of dense or empty % rows and columns ignored by COLAMD and stats (3) is the number of % garbage collections performed on the internal data structure used by % COLAMD (roughly of size 2.2*nnz(S) + 4*m + 7*n integers). % % MATLAB built-in functions are intended to generate valid sparse matrices, % with no duplicate entries, with ascending row indices of the nonzeros % in each column, with a non-negative number of entries in each column (!) % and so on. If a matrix is invalid, then COLAMD may or may not be able % to continue. If there are duplicate entries (a row index appears two or % more times in the same column) or if the row indices in a column are out % of order, then COLAMD can correct these errors by ignoring the duplicate % entries and sorting each column of its internal copy of the matrix S (the % input matrix S is not repaired, however). If a matrix is invalid in other % ways then COLAMD cannot continue, an error message is printed, and no % output arguments (P or stats) are returned. COLAMD is thus a simple way % to check a sparse matrix to see if it's valid. % % stats (4:7) provide information if COLAMD was able to continue. The % matrix is OK if stats (4) is zero, or 1 if invalid. stats (5) is the % rightmost column index that is unsorted or contains duplicate entries, % or zero if no such column exists. stats (6) is the last seen duplicate % or out-of-order row index in the column index given by stats (5), or zero % if no such row index exists. stats (7) is the number of duplicate or % out-of-order row indices. % % stats (8:20) is always zero in the current version of COLAMD (reserved % for future use). % % The ordering is followed by a column elimination tree post-ordering. % % Authors: % % The authors of the code itself are Stefan I. Larimore and Timothy A. % Davis (davis@cise.ufl.edu), University of Florida. The algorithm was % developed in collaboration with John Gilbert, Xerox PARC, and Esmond % Ng, Oak Ridge National Laboratory. % % Date: % % September 8, 2003. Version 2.3. % % Acknowledgements: % % This work was supported by the National Science Foundation, under % grants DMS-9504974 and DMS-9803599. % Notice: % % Copyright (c) 1998-2003 by the University of Florida. % All Rights Reserved. % % See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c % file) for the License. % % Availability: % % The colamd/symamd library is available at % % http://www.cise.ufl.edu/research/sparse/colamd/ % %------------------------------------------------------------------------------- % Perform the colamd ordering: %------------------------------------------------------------------------------- if (nargout <= 1 & nargin == 1) p = colamdmex (S) ; elseif (nargout <= 1 & nargin == 2) p = colamdmex (S, knobs) ; elseif (nargout == 2 & nargin == 1) [p, stats] = colamdmex (S) ; elseif (nargout == 2 & nargin == 2) [p, stats] = colamdmex (S, knobs) ; else error ('colamd: incorrect number of input and/or output arguments') ; end %------------------------------------------------------------------------------- % column elimination tree post-ordering: %------------------------------------------------------------------------------- [ignore, q] = sparsfun ('coletree', S (:,p)) ; p = p (q) ; SHAR_EOF fi # end of overwriting check if test -f 'colamdmex.c' then echo shar: will not over-write existing file "'colamdmex.c'" else cat << "SHAR_EOF" > 'colamdmex.c' /* ========================================================================== */ /* === colamd mexFunction =================================================== */ /* ========================================================================== */ /* Usage: P = colamd (A) ; P = colamd (A, knobs) ; [ P, stats ] = colamd (A) ; [ P, stats ] = colamd (A, knobs) ; Returns a permutation vector P such that the LU factorization of A (:,P) tends to be sparser than that of A. The Cholesky factorization of (A (:,P))'*(A (:,P)) will also tend to be sparser than that of A'*A. This routine provides the same functionality as COLMMD, but is much faster and returns a better permutation vector. Note that the COLMMD m-file in MATLAB 5.2 also performs a column elimination tree post-ordering. This mexFunction does not do this post-ordering. This mexFunction is a replacement for the p = sparsfun ('colmmd', A) operation. The knobs and stats vectors are optional: knobs (1) rows with more than (knobs (1))*n_col entries are removed prior to ordering. If knobs is not present, then the default is used (0.5). knobs (2) columns with more than (knobs (2))*n_row entries are removed prior to ordering, and placed last in the column permutation. If knobs is not present, then the default is used (0.5). knobs (3) print level, similar to spparms ('spumoni') stats (1) the number of dense (or empty) rows ignored stats (2) the number of dense (or empty) columms. These are ordered last, in their natural order. stats (3) the number of garbage collections performed. stats (4) return status: 0: matrix is a valid MATLAB matrix. 1: matrix has duplicate entries or unsorted columns. This should not occur in a valid MATLAB matrix, but the ordering proceeded anyway by sorting the row indices in each column and by ignoring the duplicate row indices in each column. See stats (5:7) for more information. stats (5) highest numbered column that is unsorted or has duplicate entries (zero if none) stats (6) last seen duplicate or unsorted row index (zero if none) stats (7) number of duplicate or unsorted row indices Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Notice: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c file) for the License. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/colamdmex.c file. It requires the colamd.c and colamd.h files. */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "colamd.h" #include "mex.h" #include "matrix.h" #include #include /* ========================================================================== */ /* === colamd mexFunction =================================================== */ /* ========================================================================== */ void mexFunction ( /* === Parameters ======================================================= */ int nlhs, /* number of left-hand sides */ mxArray *plhs [], /* left-hand side matrices */ int nrhs, /* number of right--hand sides */ const mxArray *prhs [] /* right-hand side matrices */ ) { /* === Local variables ================================================== */ int *A ; /* colamd's copy of the matrix, and workspace */ int *p ; /* colamd's copy of the column pointers */ int Alen ; /* size of A */ int n_col ; /* number of columns of A */ int n_row ; /* number of rows of A */ int nnz ; /* number of entries in A */ int full ; /* TRUE if input matrix full, FALSE if sparse */ double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */ double *out_perm ; /* output permutation vector */ double *out_stats ; /* output stats vector */ double *in_knobs ; /* input knobs vector */ int i ; /* loop counter */ mxArray *Ainput ; /* input matrix handle */ int spumoni ; /* verbosity variable */ int stats [COLAMD_STATS] ; /* stats for colamd */ /* === Check inputs ===================================================== */ if (nrhs < 1 || nrhs > 2 || nlhs < 0 || nlhs > 2) { mexErrMsgTxt ( "colamd: incorrect number of input and/or output arguments") ; } /* === Get knobs ======================================================== */ colamd_set_defaults (knobs) ; spumoni = 0 ; /* check for user-passed knobs */ if (nrhs == 2) { in_knobs = mxGetPr (prhs [1]) ; i = mxGetNumberOfElements (prhs [1]) ; if (i > 0) knobs [COLAMD_DENSE_ROW] = in_knobs [COLAMD_DENSE_ROW] ; if (i > 1) knobs [COLAMD_DENSE_COL] = in_knobs [COLAMD_DENSE_COL] ; if (i > 2) spumoni = (int) in_knobs [2] ; } /* print knob settings if spumoni is set */ if (spumoni > 0) { mexPrintf ("colamd: dense row fraction: %f\n", knobs [COLAMD_DENSE_ROW]) ; mexPrintf ("colamd: dense col fraction: %f\n", knobs [COLAMD_DENSE_COL]) ; } /* === If A is full, convert to a sparse matrix ========================= */ Ainput = (mxArray *) prhs [0] ; if (mxGetNumberOfDimensions (Ainput) != 2) { mexErrMsgTxt ("colamd: input matrix must be 2-dimensional") ; } full = !mxIsSparse (Ainput) ; if (full) { mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ; } /* === Allocate workspace for colamd ==================================== */ /* get size of matrix */ n_row = mxGetM (Ainput) ; n_col = mxGetN (Ainput) ; /* get column pointer vector so we can find nnz */ p = (int *) mxCalloc (n_col+1, sizeof (int)) ; (void) memcpy (p, mxGetJc (Ainput), (n_col+1)*sizeof (int)) ; nnz = p [n_col] ; Alen = colamd_recommended (nnz, n_row, n_col) ; /* === Copy input matrix into workspace ================================= */ A = (int *) mxCalloc (Alen, sizeof (int)) ; (void) memcpy (A, mxGetIr (Ainput), nnz*sizeof (int)) ; if (full) { mxDestroyArray (Ainput) ; } /* === Order the columns (destroys A) =================================== */ if (!colamd (n_row, n_col, Alen, A, p, knobs, stats)) { colamd_report (stats) ; mexErrMsgTxt ("colamd error!") ; } mxFree (A) ; /* === Return the permutation vector ==================================== */ plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; for (i = 0 ; i < n_col ; i++) { /* colamd is 0-based, but MATLAB expects this to be 1-based */ out_perm [i] = p [i] + 1 ; } mxFree (p) ; /* === Return the stats vector ========================================== */ /* print stats if spumoni > 0 */ if (spumoni > 0) { colamd_report (stats) ; } if (nlhs == 2) { plhs [1] = mxCreateDoubleMatrix (1, COLAMD_STATS, mxREAL) ; out_stats = mxGetPr (plhs [1]) ; for (i = 0 ; i < COLAMD_STATS ; i++) { out_stats [i] = stats [i] ; } /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ /* note that this correction doesn't occur if symamd returns FALSE */ out_stats [COLAMD_INFO1] ++ ; out_stats [COLAMD_INFO2] ++ ; } } SHAR_EOF fi # end of overwriting check if test -f 'symamd.m' then echo shar: will not over-write existing file "'symamd.m'" else cat << "SHAR_EOF" > 'symamd.m' function [p, stats] = symamd (S, knobs) %SYMAMD Symmetric approximate minimum degree permutation. % P = SYMAMD (S) for a symmetric positive definite matrix S, returns the % permutation vector p such that S(p,p) tends to have a sparser Cholesky % factor than S. Sometimes SYMAMD works well for symmetric indefinite % matrices too. SYMAMD tends to be faster than SYMMMD and tends to return % a better ordering. The matrix S is assumed to be symmetric; only the % strictly lower triangular part is referenced. S must be square. % % See also COLMMD, COLPERM, SPPARMS, COLAMD, SYMMMD, SYMRCM. % % Usage: P = symamd (S) % P = symamd (S, knobs) % [P, stats] = symamd (S) % [P, stats] = symamd (S, knobs) % % knobs is an optional input argument. If S is n-by-n, then rows and % columns with more than knobs(1)*n entries are removed prior to ordering, % and ordered last in the output permutation P. If the knobs parameter is not % present, then the default of 0.5 is used instead. knobs (2) controls the % printing of statistics and error messages. % % stats is an optional 20-element output vector that provides data about the % ordering and the validity of the input matrix S. Ordering statistics are % in stats (1:3). stats (1) = stats (2) is the number of dense or empty % rows and columns ignored by SYMAMD and stats (3) is the number of % garbage collections performed on the internal data structure used by % SYMAMD (roughly of size 8.4*nnz(tril(S,-1)) + 9*n integers). % % MATLAB built-in functions are intended to generate valid sparse matrices, % with no duplicate entries, with ascending row indices of the nonzeros % in each column, with a non-negative number of entries in each column (!) % and so on. If a matrix is invalid, then SYMAMD may or may not be able % to continue. If there are duplicate entries (a row index appears two or % more times in the same column) or if the row indices in a column are out % of order, then SYMAMD can correct these errors by ignoring the duplicate % entries and sorting each column of its internal copy of the matrix S (the % input matrix S is not repaired, however). If a matrix is invalid in other % ways then SYMAMD cannot continue, an error message is printed, and no % output arguments (P or stats) are returned. SYMAMD is thus a simple way % to check a sparse matrix to see if it's valid. % % stats (4:7) provide information if SYMAMD was able to continue. The % matrix is OK if stats (4) is zero, or 1 if invalid. stats (5) is the % rightmost column index that is unsorted or contains duplicate entries, % or zero if no such column exists. stats (6) is the last seen duplicate % or out-of-order row index in the column index given by stats (5), or zero % if no such row index exists. stats (7) is the number of duplicate or % out-of-order row indices. % % stats (8:20) is always zero in the current version of SYMAMD (reserved % for future use). % % The ordering is followed by a column elimination tree post-ordering. % % Authors: % % The authors of the code itself are Stefan I. Larimore and Timothy A. % Davis (davis@cise.ufl.edu), University of Florida. The algorithm was % developed in collaboration with John Gilbert, Xerox PARC, and Esmond % Ng, Oak Ridge National Laboratory. % % Date: % % September 8, 2003. Version 2.3. % % Acknowledgements: % % This work was supported by the National Science Foundation, under % grants DMS-9504974 and DMS-9803599. % Notice: % Copyright (c) 1998-2003 by the University of Florida. % All Rights Reserved. % % See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c % file) for the License. % % Availability: % % The colamd/symamd library is available at % % http://www.cise.ufl.edu/research/sparse/colamd/ % %------------------------------------------------------------------------------- % perform the symamd ordering: %------------------------------------------------------------------------------- if (nargout <= 1 & nargin == 1) p = symamdmex (S) ; elseif (nargout <= 1 & nargin == 2) p = symamdmex (S, knobs) ; elseif (nargout == 2 & nargin == 1) [p, stats] = symamdmex (S) ; elseif (nargout == 2 & nargin == 2) [p, stats] = symamdmex (S, knobs) ; else error ('symamd: incorrect number of input and/or output arguments') ; end %------------------------------------------------------------------------------- % symmetric elimination tree post-ordering: %------------------------------------------------------------------------------- [ignore, q] = sparsfun ('symetree', S (p,p)) ; p = p (q) ; SHAR_EOF fi # end of overwriting check if test -f 'symamdmex.c' then echo shar: will not over-write existing file "'symamdmex.c'" else cat << "SHAR_EOF" > 'symamdmex.c' /* ========================================================================== */ /* === symamd mexFunction =================================================== */ /* ========================================================================== */ /* Usage: P = symamd (A) ; P = symamd (A, knobs) ; [ P, stats ] = symamd (A) ; [ P, stats ] = symamd (A, knobs) ; Returns a permutation vector P such that the Cholesky factorization of A (P,P) tends to be sparser than that of A. This routine provides the same functionality as SYMMMD, but tends to be much faster and tends to return a better permutation vector. Note that the SYMMMD m-file in MATLAB 5.2 also performs a symmetric elimination tree post-ordering. This mexFunction does not do this post-ordering. This mexFunction is a replacement for the p = sparsfun ('symmmd', A) operation. A must be square, and is assummed to have a symmetric nonzero pattern. Only the nonzero pattern of the lower triangular portion of A is accessed. This routine constructs a matrix M such that the nonzero pattern of M'M is equal to A (assuming A has symmetric pattern), and then performs a column ordering of M using colamd. The knobs and stats vectors are optional: knobs (1) rows and columns with more than (knobs(1))*n entries are removed prior to ordering, and placed last in the output ordering. If knobs is not present, then the default of 0.5 is used. knobs (2) print level, similar to spparms ('spumoni') stats (1) the number of dense (or empty) rows and columms. These are ordered last, in their natural order. stats (2) (same as stats (1)) stats (3) the number of garbage collections performed. stats (4) return status: 0: matrix is a valid MATLAB matrix. 1: matrix has duplicate entries or unsorted columns. This should not occur in a valid MATLAB matrix, but the ordering proceeded anyway by ignoring the duplicate row indices in each column. See stats (5:7) for more information. stats (5) highest numbered column that is unsorted or has duplicate entries (zero if none) stats (6) last seen duplicate or unsorted row index (zero if none) stats (7) number of duplicate or unsorted row indices Authors: The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida. The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory. Date: September 8, 2003. Version 2.3. Acknowledgements: This work was supported by the National Science Foundation, under grants DMS-9504974 and DMS-9803599. Notice: Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c file) for the License. Availability: The colamd/symamd library is available at http://www.cise.ufl.edu/research/sparse/colamd/ This is the http://www.cise.ufl.edu/research/sparse/colamd/symamdmex.c file. It requires the colamd.c and colamd.h files. */ /* ========================================================================== */ /* === Include files ======================================================== */ /* ========================================================================== */ #include "colamd.h" #include "mex.h" #include "matrix.h" #include /* ========================================================================== */ /* === symamd mexFunction =================================================== */ /* ========================================================================== */ void mexFunction ( /* === Parameters ======================================================= */ int nlhs, /* number of left-hand sides */ mxArray *plhs [], /* left-hand side matrices */ int nrhs, /* number of right--hand sides */ const mxArray *prhs [] /* right-hand side matrices */ ) { /* === Local variables ================================================== */ int *perm ; /* column ordering of M and ordering of A */ int *A ; /* row indices of input matrix A */ int *p ; /* column pointers of input matrix A */ int n_col ; /* number of columns of A */ int n_row ; /* number of rows of A */ int full ; /* TRUE if input matrix full, FALSE if sparse */ double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */ double *out_perm ; /* output permutation vector */ double *out_stats ; /* output stats vector */ double *in_knobs ; /* input knobs vector */ int i ; /* loop counter */ mxArray *Ainput ; /* input matrix handle */ int spumoni ; /* verbosity variable */ int stats [COLAMD_STATS] ; /* stats for symamd */ /* === Check inputs ===================================================== */ if (nrhs < 1 || nrhs > 2 || nlhs < 0 || nlhs > 2) { mexErrMsgTxt ( "symamd: incorrect number of input and/or output arguments.") ; } /* === Get knobs ======================================================== */ colamd_set_defaults (knobs) ; spumoni = 0 ; /* check for user-passed knobs */ if (nrhs == 2) { in_knobs = mxGetPr (prhs [1]) ; i = mxGetNumberOfElements (prhs [1]) ; if (i > 0) knobs [COLAMD_DENSE_ROW] = in_knobs [COLAMD_DENSE_ROW] ; if (i > 1) spumoni = (int) in_knobs [1] ; } /* print knob settings if spumoni > 0 */ if (spumoni > 0) { mexPrintf ("symamd: dense row/col fraction: %f\n", knobs [COLAMD_DENSE_ROW]) ; } /* === If A is full, convert to a sparse matrix ========================= */ Ainput = (mxArray *) prhs [0] ; if (mxGetNumberOfDimensions (Ainput) != 2) { mexErrMsgTxt ("symamd: input matrix must be 2-dimensional.") ; } full = !mxIsSparse (Ainput) ; if (full) { mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ; } /* === Allocate workspace for symamd ==================================== */ /* get size of matrix */ n_row = mxGetM (Ainput) ; n_col = mxGetN (Ainput) ; if (n_col != n_row) { mexErrMsgTxt ("symamd: matrix must be square.") ; } A = mxGetIr (Ainput) ; p = mxGetJc (Ainput) ; perm = (int *) mxCalloc (n_col+1, sizeof (int)) ; /* === Order the rows and columns of A (does not destroy A) ============= */ if (!symamd (n_col, A, p, perm, knobs, stats, &mxCalloc, &mxFree)) { symamd_report (stats) ; mexErrMsgTxt ("symamd error!") ; } if (full) { mxDestroyArray (Ainput) ; } /* === Return the permutation vector ==================================== */ plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ; out_perm = mxGetPr (plhs [0]) ; for (i = 0 ; i < n_col ; i++) { /* symamd is 0-based, but MATLAB expects this to be 1-based */ out_perm [i] = perm [i] + 1 ; } mxFree (perm) ; /* === Return the stats vector ========================================== */ /* print stats if spumoni > 0 */ if (spumoni > 0) { symamd_report (stats) ; } if (nlhs == 2) { plhs [1] = mxCreateDoubleMatrix (1, COLAMD_STATS, mxREAL) ; out_stats = mxGetPr (plhs [1]) ; for (i = 0 ; i < COLAMD_STATS ; i++) { out_stats [i] = stats [i] ; } /* fix stats (5) and (6), for 1-based information on jumbled matrix. */ /* note that this correction doesn't occur if symamd returns FALSE */ out_stats [COLAMD_INFO1] ++ ; out_stats [COLAMD_INFO2] ++ ; } } SHAR_EOF fi # end of overwriting check cd .. cd .. if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << "SHAR_EOF" > 'README' The COLAMD ordering method - Version 2.3 ------------------------------------------------------------------------------- The COLAMD column approximate minimum degree ordering algorithm computes a permutation vector P such that the LU factorization of A (:,P) tends to be sparser than that of A. The Cholesky factorization of (A (:,P))'*(A (:,P)) will also tend to be sparser than that of A'*A. SYMAMD is a symmetric minimum degree ordering method based on COLAMD, available as a MATLAB-callable function. It constructs a matrix M such that M'*M has the same pattern as A, and then uses COLAMD to compute a column ordering of M. Colamd and symamd tend to be faster and generate better orderings than their MATLAB counterparts, colmmd and symmmd. To compile and test the colamd m-files and mexFunctions, just unpack the colamd2.3/ directory from the colamd2.3.tar.gz file, and run MATLAB from within that directory. Next, type colamd_test to compile and test colamd and symamd. This will work on any computer with MATLAB (Unix, PC, or Mac). Alternatively, type "make" (in Unix) to compile and run a simple example C code, without using MATLAB. Colamd 2.0 is a built-in routine in MATLAB V6.0, available from The Mathworks, Inc. Under most cases, the compiled codes from Versions 2.0 through 2.2 do not differ. Colamd Versions 2.2 and 2.3 differ only in their mexFunction interaces to MATLAB. To use colamd and symamd within an application written in C, all you need are colamd.c and colamd.h, which are the C-callable colamd/symamd codes. See colamd.c for more information on how to call colamd from a C program. Copyright (c) 1998-2003 by the University of Florida. All Rights Reserved. See http://www.cise.ufl.edu/research/sparse/colamd (the colamd.c file) for the License. Related papers: "A column approximate minimum degree ordering algorithm", Timothy A. Davis, John R. Gilbert, Stefan I. Larimore, and Esmond G. Ng. ACM Trans. on Mathematical Software. "Algorithm 8xx: COLAMD, a column approximate minimum degree ordering algorithm", Timothy A. Davis, John R. Gilbert, Stefan I. Larimore, and Esmond G. Ng. ACM Trans. on Mathematical Software. "An approximate minimum degree column ordering algorithm", S. I. Larimore, MS Thesis, Dept. of Computer and Information Science and Engineering, University of Florida, Gainesville, FL, 1998. CISE Tech Report TR-98-016. Available at ftp://ftp.cise.ufl.edu/cis/tech-reports/tr98/tr98-016.ps via anonymous ftp. Approximate Deficiency for Ordering the Columns of a Matrix, J. L. Kern, Senior Thesis, Dept. of Computer and Information Science and Engineering, University of Florida, Gainesville, FL, 1999. Available at http://www.cise.ufl.edu/~davis/Kern/kern.ps Authors: Stefan I. Larimore and Timothy A. Davis, University of Florida, in collaboration with John Gilbert, Xerox PARC (now at UC Santa Barbara), and Esmong Ng, Lawrence Berkeley National Laboratory (much of this work he did while at Oak Ridge National Laboratory). COLAMD files (Version 2.3, September 8, 2003): colamd2.3.tar.gz: All files, as a gzipped, Unix tar file. The *.m, and *mex.c files are for use in MATLAB. colamd.c: the primary colamd computational kernel. colamd.h: include file for colamd/symamd library. colamd.m: the MATLAB interface to colamd. colamd_demo.m: MATLAB demo file for colamd and symamd (also compiles the colamdmex and symamdmex mexFunctions). colamdmex.c: colamd mexFunction for use in MATLAB. colamd_example.c: example C main program that calls colamd and symamd. colamd_example.out: output of colamd_example.c. Makefile: Makefile for colamd_example.c symamd.m: the MATLAB interface to symamd. symamdmex.c: symamd mexFunction for use in MATLAB. README: this file ChangeLog: a log of changes since Version 1.0. colamd_test.m: test code colamdtestmex.c: test code luflops.m: test code symamdtestmex.c: test code SHAR_EOF fi # end of overwriting check # End of shell archive exit 0