alpar@9: /* ========================================================================== */ alpar@9: /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* COLAMD / SYMAMD alpar@9: alpar@9: colamd: an approximate minimum degree column ordering algorithm, alpar@9: for LU factorization of symmetric or unsymmetric matrices, alpar@9: QR factorization, least squares, interior point methods for alpar@9: linear programming problems, and other related problems. alpar@9: alpar@9: symamd: an approximate minimum degree ordering algorithm for Cholesky alpar@9: factorization of symmetric matrices. alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Colamd computes a permutation Q such that the Cholesky factorization of alpar@9: (AQ)'(AQ) has less fill-in and requires fewer floating point operations alpar@9: than A'A. This also provides a good ordering for sparse partial alpar@9: pivoting methods, P(AQ) = LU, where Q is computed prior to numerical alpar@9: factorization, and P is computed during numerical factorization via alpar@9: conventional partial pivoting with row interchanges. Colamd is the alpar@9: column ordering method used in SuperLU, part of the ScaLAPACK library. alpar@9: It is also available as built-in function in MATLAB Version 6, alpar@9: available from MathWorks, Inc. (http://www.mathworks.com). This alpar@9: routine can be used in place of colmmd in MATLAB. alpar@9: alpar@9: Symamd computes a permutation P of a symmetric matrix A such that the alpar@9: Cholesky factorization of PAP' has less fill-in and requires fewer alpar@9: floating point operations than A. Symamd constructs a matrix M such alpar@9: that M'M has the same nonzero pattern of A, and then orders the columns alpar@9: of M using colmmd. The column ordering of M is then returned as the alpar@9: row and column ordering P of A. alpar@9: alpar@9: Authors: alpar@9: alpar@9: The authors of the code itself are Stefan I. Larimore and Timothy A. alpar@9: Davis (davis at cise.ufl.edu), University of Florida. The algorithm was alpar@9: developed in collaboration with John Gilbert, Xerox PARC, and Esmond alpar@9: Ng, Oak Ridge National Laboratory. alpar@9: alpar@9: Acknowledgements: alpar@9: alpar@9: This work was supported by the National Science Foundation, under alpar@9: grants DMS-9504974 and DMS-9803599. alpar@9: alpar@9: Copyright and License: alpar@9: alpar@9: Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved. alpar@9: COLAMD is also available under alternate licenses, contact T. Davis alpar@9: for details. alpar@9: alpar@9: This library is free software; you can redistribute it and/or alpar@9: modify it under the terms of the GNU Lesser General Public alpar@9: License as published by the Free Software Foundation; either alpar@9: version 2.1 of the License, or (at your option) any later version. alpar@9: alpar@9: This library is distributed in the hope that it will be useful, alpar@9: but WITHOUT ANY WARRANTY; without even the implied warranty of alpar@9: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU alpar@9: Lesser General Public License for more details. alpar@9: alpar@9: You should have received a copy of the GNU Lesser General Public alpar@9: License along with this library; if not, write to the Free Software alpar@9: Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 alpar@9: USA alpar@9: alpar@9: Permission is hereby granted to use or copy this program under the alpar@9: terms of the GNU LGPL, provided that the Copyright, this License, alpar@9: and the Availability of the original version is retained on all copies. alpar@9: User documentation of any code that uses this code or any modified alpar@9: version of this code must cite the Copyright, this License, the alpar@9: Availability note, and "Used by permission." Permission to modify alpar@9: the code and to distribute modified code is granted, provided the alpar@9: Copyright, this License, and the Availability note are retained, alpar@9: and a notice that the code was modified is included. alpar@9: alpar@9: Availability: alpar@9: alpar@9: The colamd/symamd library is available at alpar@9: alpar@9: http://www.cise.ufl.edu/research/sparse/colamd/ alpar@9: alpar@9: This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c alpar@9: file. It requires the colamd.h file. It is required by the colamdmex.c alpar@9: and symamdmex.c files, for the MATLAB interface to colamd and symamd. alpar@9: Appears as ACM Algorithm 836. alpar@9: alpar@9: See the ChangeLog file for changes since Version 1.0. alpar@9: alpar@9: References: alpar@9: alpar@9: T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column alpar@9: minimum degree ordering algorithm, ACM Transactions on Mathematical alpar@9: Software, vol. 30, no. 3., pp. 353-376, 2004. alpar@9: alpar@9: T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD, alpar@9: an approximate column minimum degree ordering algorithm, ACM alpar@9: Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380, alpar@9: 2004. alpar@9: alpar@9: */ alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Description of user-callable routines ================================ */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* COLAMD includes both int and UF_long versions of all its routines. The alpar@9: * description below is for the int version. For UF_long, all int arguments alpar@9: * become UF_long. UF_long is normally defined as long, except for WIN64. alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: colamd_recommended: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: size_t colamd_recommended (int nnz, int n_row, int n_col) ; alpar@9: size_t colamd_l_recommended (UF_long nnz, UF_long n_row, alpar@9: UF_long n_col) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Returns recommended value of Alen for use by colamd. Returns 0 alpar@9: if any input argument is negative. The use of this routine alpar@9: is optional. Not needed for symamd, which dynamically allocates alpar@9: its own memory. alpar@9: alpar@9: Note that in v2.4 and earlier, these routines returned int or long. alpar@9: They now return a value of type size_t. alpar@9: alpar@9: Arguments (all input arguments): alpar@9: alpar@9: int nnz ; Number of nonzeros in the matrix A. This must alpar@9: be the same value as p [n_col] in the call to alpar@9: colamd - otherwise you will get a wrong value alpar@9: of the recommended memory to use. alpar@9: alpar@9: int n_row ; Number of rows in the matrix A. alpar@9: alpar@9: int n_col ; Number of columns in the matrix A. alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: colamd_set_defaults: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: colamd_set_defaults (double knobs [COLAMD_KNOBS]) ; alpar@9: colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Sets the default parameters. The use of this routine is optional. alpar@9: alpar@9: Arguments: alpar@9: alpar@9: double knobs [COLAMD_KNOBS] ; Output only. alpar@9: alpar@9: NOTE: the meaning of the dense row/col knobs has changed in v2.4 alpar@9: alpar@9: knobs [0] and knobs [1] control dense row and col detection: alpar@9: alpar@9: Colamd: rows with more than alpar@9: max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col)) alpar@9: entries are removed prior to ordering. Columns with more than alpar@9: max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col))) alpar@9: entries are removed prior to alpar@9: ordering, and placed last in the output column ordering. alpar@9: alpar@9: Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0]. alpar@9: Rows and columns with more than alpar@9: max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n)) alpar@9: entries are removed prior to ordering, and placed last in the alpar@9: output ordering. alpar@9: alpar@9: COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1, alpar@9: respectively, in colamd.h. Default values of these two knobs alpar@9: are both 10. Currently, only knobs [0] and knobs [1] are alpar@9: used, but future versions may use more knobs. If so, they will alpar@9: be properly set to their defaults by the future version of alpar@9: colamd_set_defaults, so that the code that calls colamd will alpar@9: not need to change, assuming that you either use alpar@9: colamd_set_defaults, or pass a (double *) NULL pointer as the alpar@9: knobs array to colamd or symamd. alpar@9: alpar@9: knobs [2]: aggressive absorption alpar@9: alpar@9: knobs [COLAMD_AGGRESSIVE] controls whether or not to do alpar@9: aggressive absorption during the ordering. Default is TRUE. alpar@9: alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: colamd: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: int colamd (int n_row, int n_col, int Alen, int *A, int *p, alpar@9: double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ; alpar@9: UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen, alpar@9: UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS], alpar@9: UF_long stats [COLAMD_STATS]) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Computes a column ordering (Q) of A such that P(AQ)=LU or alpar@9: (AQ)'AQ=LL' have less fill-in and require fewer floating point alpar@9: operations than factorizing the unpermuted matrix A or A'A, alpar@9: respectively. alpar@9: alpar@9: Returns: alpar@9: alpar@9: TRUE (1) if successful, FALSE (0) otherwise. alpar@9: alpar@9: Arguments: alpar@9: alpar@9: int n_row ; Input argument. alpar@9: alpar@9: Number of rows in the matrix A. alpar@9: Restriction: n_row >= 0. alpar@9: Colamd returns FALSE if n_row is negative. alpar@9: alpar@9: int n_col ; Input argument. alpar@9: alpar@9: Number of columns in the matrix A. alpar@9: Restriction: n_col >= 0. alpar@9: Colamd returns FALSE if n_col is negative. alpar@9: alpar@9: int Alen ; Input argument. alpar@9: alpar@9: Restriction (see note): alpar@9: Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col alpar@9: Colamd returns FALSE if these conditions are not met. alpar@9: alpar@9: Note: this restriction makes an modest assumption regarding alpar@9: the size of the two typedef's structures in colamd.h. alpar@9: We do, however, guarantee that alpar@9: alpar@9: Alen >= colamd_recommended (nnz, n_row, n_col) alpar@9: alpar@9: will be sufficient. Note: the macro version does not check alpar@9: for integer overflow, and thus is not recommended. Use alpar@9: the colamd_recommended routine instead. alpar@9: alpar@9: int A [Alen] ; Input argument, undefined on output. alpar@9: alpar@9: A is an integer array of size Alen. Alen must be at least as alpar@9: large as the bare minimum value given above, but this is very alpar@9: low, and can result in excessive run time. For best alpar@9: performance, we recommend that Alen be greater than or equal to alpar@9: colamd_recommended (nnz, n_row, n_col), which adds alpar@9: nnz/5 to the bare minimum value given above. alpar@9: alpar@9: On input, the row indices of the entries in column c of the alpar@9: matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices alpar@9: in a given column c need not be in ascending order, and alpar@9: duplicate row indices may be be present. However, colamd will alpar@9: work a little faster if both of these conditions are met alpar@9: (Colamd puts the matrix into this format, if it finds that the alpar@9: the conditions are not met). alpar@9: alpar@9: The matrix is 0-based. That is, rows are in the range 0 to alpar@9: n_row-1, and columns are in the range 0 to n_col-1. Colamd alpar@9: returns FALSE if any row index is out of range. alpar@9: alpar@9: The contents of A are modified during ordering, and are alpar@9: undefined on output. alpar@9: alpar@9: int p [n_col+1] ; Both input and output argument. alpar@9: alpar@9: p is an integer array of size n_col+1. On input, it holds the alpar@9: "pointers" for the column form of the matrix A. Column c of alpar@9: the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first alpar@9: entry, p [0], must be zero, and p [c] <= p [c+1] must hold alpar@9: for all c in the range 0 to n_col-1. The value p [n_col] is alpar@9: thus the total number of entries in the pattern of the matrix A. alpar@9: Colamd returns FALSE if these conditions are not met. alpar@9: alpar@9: On output, if colamd returns TRUE, the array p holds the column alpar@9: permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is alpar@9: the first column index in the new ordering, and p [n_col-1] is alpar@9: the last. That is, p [k] = j means that column j of A is the alpar@9: kth pivot column, in AQ, where k is in the range 0 to n_col-1 alpar@9: (p [0] = j means that column j of A is the first column in AQ). alpar@9: alpar@9: If colamd returns FALSE, then no permutation is returned, and alpar@9: p is undefined on output. alpar@9: alpar@9: double knobs [COLAMD_KNOBS] ; Input argument. alpar@9: alpar@9: See colamd_set_defaults for a description. alpar@9: alpar@9: int stats [COLAMD_STATS] ; Output argument. alpar@9: alpar@9: Statistics on the ordering, and error status. alpar@9: See colamd.h for related definitions. alpar@9: Colamd returns FALSE if stats is not present. alpar@9: alpar@9: stats [0]: number of dense or empty rows ignored. alpar@9: alpar@9: stats [1]: number of dense or empty columns ignored (and alpar@9: ordered last in the output permutation p) alpar@9: Note that a row can become "empty" if it alpar@9: contains only "dense" and/or "empty" columns, alpar@9: and similarly a column can become "empty" if it alpar@9: only contains "dense" and/or "empty" rows. alpar@9: alpar@9: stats [2]: number of garbage collections performed. alpar@9: This can be excessively high if Alen is close alpar@9: to the minimum required value. alpar@9: alpar@9: stats [3]: status code. < 0 is an error code. alpar@9: > 1 is a warning or notice. alpar@9: alpar@9: 0 OK. Each column of the input matrix contained alpar@9: row indices in increasing order, with no alpar@9: duplicates. alpar@9: alpar@9: 1 OK, but columns of input matrix were jumbled alpar@9: (unsorted columns or duplicate entries). Colamd alpar@9: had to do some extra work to sort the matrix alpar@9: first and remove duplicate entries, but it alpar@9: still was able to return a valid permutation alpar@9: (return value of colamd was TRUE). alpar@9: alpar@9: stats [4]: highest numbered column that alpar@9: is unsorted or has duplicate alpar@9: entries. alpar@9: stats [5]: last seen duplicate or alpar@9: unsorted row index. alpar@9: stats [6]: number of duplicate or alpar@9: unsorted row indices. alpar@9: alpar@9: -1 A is a null pointer alpar@9: alpar@9: -2 p is a null pointer alpar@9: alpar@9: -3 n_row is negative alpar@9: alpar@9: stats [4]: n_row alpar@9: alpar@9: -4 n_col is negative alpar@9: alpar@9: stats [4]: n_col alpar@9: alpar@9: -5 number of nonzeros in matrix is negative alpar@9: alpar@9: stats [4]: number of nonzeros, p [n_col] alpar@9: alpar@9: -6 p [0] is nonzero alpar@9: alpar@9: stats [4]: p [0] alpar@9: alpar@9: -7 A is too small alpar@9: alpar@9: stats [4]: required size alpar@9: stats [5]: actual size (Alen) alpar@9: alpar@9: -8 a column has a negative number of entries alpar@9: alpar@9: stats [4]: column with < 0 entries alpar@9: stats [5]: number of entries in col alpar@9: alpar@9: -9 a row index is out of bounds alpar@9: alpar@9: stats [4]: column with bad row index alpar@9: stats [5]: bad row index alpar@9: stats [6]: n_row, # of rows of matrx alpar@9: alpar@9: -10 (unused; see symamd.c) alpar@9: alpar@9: -999 (unused; see symamd.c) alpar@9: alpar@9: Future versions may return more statistics in the stats array. alpar@9: alpar@9: Example: alpar@9: alpar@9: See http://www.cise.ufl.edu/research/sparse/colamd/example.c alpar@9: for a complete example. alpar@9: alpar@9: To order the columns of a 5-by-4 matrix with 11 nonzero entries in alpar@9: the following nonzero pattern alpar@9: alpar@9: x 0 x 0 alpar@9: x 0 x x alpar@9: 0 x x 0 alpar@9: 0 0 x x alpar@9: x x 0 0 alpar@9: alpar@9: with default knobs and no output statistics, do the following: alpar@9: alpar@9: #include "colamd.h" alpar@9: #define ALEN 100 alpar@9: int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ; alpar@9: int p [ ] = {0, 3, 5, 9, 11} ; alpar@9: int stats [COLAMD_STATS] ; alpar@9: colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ; alpar@9: alpar@9: The permutation is returned in the array p, and A is destroyed. alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: symamd: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: int symamd (int n, int *A, int *p, int *perm, alpar@9: double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS], alpar@9: void (*allocate) (size_t, size_t), void (*release) (void *)) ; alpar@9: UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm, alpar@9: double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS], alpar@9: void (*allocate) (size_t, size_t), void (*release) (void *)) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: The symamd routine computes an ordering P of a symmetric sparse alpar@9: matrix A such that the Cholesky factorization PAP' = LL' remains alpar@9: sparse. It is based on a column ordering of a matrix M constructed alpar@9: so that the nonzero pattern of M'M is the same as A. The matrix A alpar@9: is assumed to be symmetric; only the strictly lower triangular part alpar@9: is accessed. You must pass your selected memory allocator (usually alpar@9: calloc/free or mxCalloc/mxFree) to symamd, for it to allocate alpar@9: memory for the temporary matrix M. alpar@9: alpar@9: Returns: alpar@9: alpar@9: TRUE (1) if successful, FALSE (0) otherwise. alpar@9: alpar@9: Arguments: alpar@9: alpar@9: int n ; Input argument. alpar@9: alpar@9: Number of rows and columns in the symmetrix matrix A. alpar@9: Restriction: n >= 0. alpar@9: Symamd returns FALSE if n is negative. alpar@9: alpar@9: int A [nnz] ; Input argument. alpar@9: alpar@9: A is an integer array of size nnz, where nnz = p [n]. alpar@9: alpar@9: The row indices of the entries in column c of the matrix are alpar@9: held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a alpar@9: given column c need not be in ascending order, and duplicate alpar@9: row indices may be present. However, symamd will run faster alpar@9: if the columns are in sorted order with no duplicate entries. alpar@9: alpar@9: The matrix is 0-based. That is, rows are in the range 0 to alpar@9: n-1, and columns are in the range 0 to n-1. Symamd alpar@9: returns FALSE if any row index is out of range. alpar@9: alpar@9: The contents of A are not modified. alpar@9: alpar@9: int p [n+1] ; Input argument. alpar@9: alpar@9: p is an integer array of size n+1. On input, it holds the alpar@9: "pointers" for the column form of the matrix A. Column c of alpar@9: the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first alpar@9: entry, p [0], must be zero, and p [c] <= p [c+1] must hold alpar@9: for all c in the range 0 to n-1. The value p [n] is alpar@9: thus the total number of entries in the pattern of the matrix A. alpar@9: Symamd returns FALSE if these conditions are not met. alpar@9: alpar@9: The contents of p are not modified. alpar@9: alpar@9: int perm [n+1] ; Output argument. alpar@9: alpar@9: On output, if symamd returns TRUE, the array perm holds the alpar@9: permutation P, where perm [0] is the first index in the new alpar@9: ordering, and perm [n-1] is the last. That is, perm [k] = j alpar@9: means that row and column j of A is the kth column in PAP', alpar@9: where k is in the range 0 to n-1 (perm [0] = j means alpar@9: that row and column j of A are the first row and column in alpar@9: PAP'). The array is used as a workspace during the ordering, alpar@9: which is why it must be of length n+1, not just n. alpar@9: alpar@9: double knobs [COLAMD_KNOBS] ; Input argument. alpar@9: alpar@9: See colamd_set_defaults for a description. alpar@9: alpar@9: int stats [COLAMD_STATS] ; Output argument. alpar@9: alpar@9: Statistics on the ordering, and error status. alpar@9: See colamd.h for related definitions. alpar@9: Symamd returns FALSE if stats is not present. alpar@9: alpar@9: stats [0]: number of dense or empty row and columns ignored alpar@9: (and ordered last in the output permutation alpar@9: perm). Note that a row/column can become alpar@9: "empty" if it contains only "dense" and/or alpar@9: "empty" columns/rows. alpar@9: alpar@9: stats [1]: (same as stats [0]) alpar@9: alpar@9: stats [2]: number of garbage collections performed. alpar@9: alpar@9: stats [3]: status code. < 0 is an error code. alpar@9: > 1 is a warning or notice. alpar@9: alpar@9: 0 OK. Each column of the input matrix contained alpar@9: row indices in increasing order, with no alpar@9: duplicates. alpar@9: alpar@9: 1 OK, but columns of input matrix were jumbled alpar@9: (unsorted columns or duplicate entries). Symamd alpar@9: had to do some extra work to sort the matrix alpar@9: first and remove duplicate entries, but it alpar@9: still was able to return a valid permutation alpar@9: (return value of symamd was TRUE). alpar@9: alpar@9: stats [4]: highest numbered column that alpar@9: is unsorted or has duplicate alpar@9: entries. alpar@9: stats [5]: last seen duplicate or alpar@9: unsorted row index. alpar@9: stats [6]: number of duplicate or alpar@9: unsorted row indices. alpar@9: alpar@9: -1 A is a null pointer alpar@9: alpar@9: -2 p is a null pointer alpar@9: alpar@9: -3 (unused, see colamd.c) alpar@9: alpar@9: -4 n is negative alpar@9: alpar@9: stats [4]: n alpar@9: alpar@9: -5 number of nonzeros in matrix is negative alpar@9: alpar@9: stats [4]: # of nonzeros (p [n]). alpar@9: alpar@9: -6 p [0] is nonzero alpar@9: alpar@9: stats [4]: p [0] alpar@9: alpar@9: -7 (unused) alpar@9: alpar@9: -8 a column has a negative number of entries alpar@9: alpar@9: stats [4]: column with < 0 entries alpar@9: stats [5]: number of entries in col alpar@9: alpar@9: -9 a row index is out of bounds alpar@9: alpar@9: stats [4]: column with bad row index alpar@9: stats [5]: bad row index alpar@9: stats [6]: n_row, # of rows of matrx alpar@9: alpar@9: -10 out of memory (unable to allocate temporary alpar@9: workspace for M or count arrays using the alpar@9: "allocate" routine passed into symamd). alpar@9: alpar@9: Future versions may return more statistics in the stats array. alpar@9: alpar@9: void * (*allocate) (size_t, size_t) alpar@9: alpar@9: A pointer to a function providing memory allocation. The alpar@9: allocated memory must be returned initialized to zero. For a alpar@9: C application, this argument should normally be a pointer to alpar@9: calloc. For a MATLAB mexFunction, the routine mxCalloc is alpar@9: passed instead. alpar@9: alpar@9: void (*release) (size_t, size_t) alpar@9: alpar@9: A pointer to a function that frees memory allocated by the alpar@9: memory allocation routine above. For a C application, this alpar@9: argument should normally be a pointer to free. For a MATLAB alpar@9: mexFunction, the routine mxFree is passed instead. alpar@9: alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: colamd_report: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: colamd_report (int stats [COLAMD_STATS]) ; alpar@9: colamd_l_report (UF_long stats [COLAMD_STATS]) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Prints the error status and statistics recorded in the stats alpar@9: array on the standard error output (for a standard C routine) alpar@9: or on the MATLAB output (for a mexFunction). alpar@9: alpar@9: Arguments: alpar@9: alpar@9: int stats [COLAMD_STATS] ; Input only. Statistics from colamd. alpar@9: alpar@9: alpar@9: ---------------------------------------------------------------------------- alpar@9: symamd_report: alpar@9: ---------------------------------------------------------------------------- alpar@9: alpar@9: C syntax: alpar@9: alpar@9: #include "colamd.h" alpar@9: symamd_report (int stats [COLAMD_STATS]) ; alpar@9: symamd_l_report (UF_long stats [COLAMD_STATS]) ; alpar@9: alpar@9: Purpose: alpar@9: alpar@9: Prints the error status and statistics recorded in the stats alpar@9: array on the standard error output (for a standard C routine) alpar@9: or on the MATLAB output (for a mexFunction). alpar@9: alpar@9: Arguments: alpar@9: alpar@9: int stats [COLAMD_STATS] ; Input only. Statistics from symamd. alpar@9: alpar@9: alpar@9: */ alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Scaffolding code definitions ======================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* Ensure that debugging is turned off: */ alpar@9: #ifndef NDEBUG alpar@9: #define NDEBUG alpar@9: #endif alpar@9: alpar@9: /* turn on debugging by uncommenting the following line alpar@9: #undef NDEBUG alpar@9: */ alpar@9: alpar@9: /* alpar@9: Our "scaffolding code" philosophy: In our opinion, well-written library alpar@9: code should keep its "debugging" code, and just normally have it turned off alpar@9: by the compiler so as not to interfere with performance. This serves alpar@9: several purposes: alpar@9: alpar@9: (1) assertions act as comments to the reader, telling you what the code alpar@9: expects at that point. All assertions will always be true (unless alpar@9: there really is a bug, of course). alpar@9: alpar@9: (2) leaving in the scaffolding code assists anyone who would like to modify alpar@9: the code, or understand the algorithm (by reading the debugging output, alpar@9: one can get a glimpse into what the code is doing). alpar@9: alpar@9: (3) (gasp!) for actually finding bugs. This code has been heavily tested alpar@9: and "should" be fully functional and bug-free ... but you never know... alpar@9: alpar@9: The code will become outrageously slow when debugging is alpar@9: enabled. To control the level of debugging output, set an environment alpar@9: variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging, alpar@9: you should see the following message on the standard output: alpar@9: alpar@9: colamd: debug version, D = 1 (THIS WILL BE SLOW!) alpar@9: alpar@9: or a similar message for symamd. If you don't, then debugging has not alpar@9: been enabled. alpar@9: alpar@9: */ alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Include files ======================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: #include "colamd.h" alpar@9: alpar@9: #if 0 /* by mao */ alpar@9: #include alpar@9: #include alpar@9: alpar@9: #ifdef MATLAB_MEX_FILE alpar@9: #include "mex.h" alpar@9: #include "matrix.h" alpar@9: #endif /* MATLAB_MEX_FILE */ alpar@9: alpar@9: #if !defined (NPRINT) || !defined (NDEBUG) alpar@9: #include alpar@9: #endif alpar@9: alpar@9: #ifndef NULL alpar@9: #define NULL ((void *) 0) alpar@9: #endif alpar@9: #endif alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === int or UF_long ======================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: #if 0 /* by mao */ alpar@9: /* define UF_long */ alpar@9: #include "UFconfig.h" alpar@9: #endif alpar@9: alpar@9: #ifdef DLONG alpar@9: alpar@9: #define Int UF_long alpar@9: #define ID UF_long_id alpar@9: #define Int_MAX UF_long_max alpar@9: alpar@9: #define COLAMD_recommended colamd_l_recommended alpar@9: #define COLAMD_set_defaults colamd_l_set_defaults alpar@9: #define COLAMD_MAIN colamd_l alpar@9: #define SYMAMD_MAIN symamd_l alpar@9: #define COLAMD_report colamd_l_report alpar@9: #define SYMAMD_report symamd_l_report alpar@9: alpar@9: #else alpar@9: alpar@9: #define Int int alpar@9: #define ID "%d" alpar@9: #define Int_MAX INT_MAX alpar@9: alpar@9: #define COLAMD_recommended colamd_recommended alpar@9: #define COLAMD_set_defaults colamd_set_defaults alpar@9: #define COLAMD_MAIN colamd alpar@9: #define SYMAMD_MAIN symamd alpar@9: #define COLAMD_report colamd_report alpar@9: #define SYMAMD_report symamd_report alpar@9: alpar@9: #endif alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Row and Column structures ============================================ */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* User code that makes use of the colamd/symamd routines need not directly */ alpar@9: /* reference these structures. They are used only for colamd_recommended. */ alpar@9: alpar@9: typedef struct Colamd_Col_struct alpar@9: { alpar@9: Int start ; /* index for A of first row in this column, or DEAD */ alpar@9: /* if column is dead */ alpar@9: Int length ; /* number of rows in this column */ alpar@9: union alpar@9: { alpar@9: Int thickness ; /* number of original columns represented by this */ alpar@9: /* col, if the column is alive */ alpar@9: Int parent ; /* parent in parent tree super-column structure, if */ alpar@9: /* the column is dead */ alpar@9: } shared1 ; alpar@9: union alpar@9: { alpar@9: Int score ; /* the score used to maintain heap, if col is alive */ alpar@9: Int order ; /* pivot ordering of this column, if col is dead */ alpar@9: } shared2 ; alpar@9: union alpar@9: { alpar@9: Int headhash ; /* head of a hash bucket, if col is at the head of */ alpar@9: /* a degree list */ alpar@9: Int hash ; /* hash value, if col is not in a degree list */ alpar@9: Int prev ; /* previous column in degree list, if col is in a */ alpar@9: /* degree list (but not at the head of a degree list) */ alpar@9: } shared3 ; alpar@9: union alpar@9: { alpar@9: Int degree_next ; /* next column, if col is in a degree list */ alpar@9: Int hash_next ; /* next column, if col is in a hash list */ alpar@9: } shared4 ; alpar@9: alpar@9: } Colamd_Col ; alpar@9: alpar@9: typedef struct Colamd_Row_struct alpar@9: { alpar@9: Int start ; /* index for A of first col in this row */ alpar@9: Int length ; /* number of principal columns in this row */ alpar@9: union alpar@9: { alpar@9: Int degree ; /* number of principal & non-principal columns in row */ alpar@9: Int p ; /* used as a row pointer in init_rows_cols () */ alpar@9: } shared1 ; alpar@9: union alpar@9: { alpar@9: Int mark ; /* for computing set differences and marking dead rows*/ alpar@9: Int first_column ;/* first column in row (used in garbage collection) */ alpar@9: } shared2 ; alpar@9: alpar@9: } Colamd_Row ; alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Definitions ========================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */ alpar@9: #define PUBLIC alpar@9: #define PRIVATE static alpar@9: alpar@9: #define DENSE_DEGREE(alpha,n) \ alpar@9: ((Int) MAX (16.0, (alpha) * sqrt ((double) (n)))) alpar@9: alpar@9: #define MAX(a,b) (((a) > (b)) ? (a) : (b)) alpar@9: #define MIN(a,b) (((a) < (b)) ? (a) : (b)) alpar@9: alpar@9: #define ONES_COMPLEMENT(r) (-(r)-1) alpar@9: alpar@9: /* -------------------------------------------------------------------------- */ alpar@9: /* Change for version 2.1: define TRUE and FALSE only if not yet defined */ alpar@9: /* -------------------------------------------------------------------------- */ alpar@9: alpar@9: #ifndef TRUE alpar@9: #define TRUE (1) alpar@9: #endif alpar@9: alpar@9: #ifndef FALSE alpar@9: #define FALSE (0) alpar@9: #endif alpar@9: alpar@9: /* -------------------------------------------------------------------------- */ alpar@9: alpar@9: #define EMPTY (-1) alpar@9: alpar@9: /* Row and column status */ alpar@9: #define ALIVE (0) alpar@9: #define DEAD (-1) alpar@9: alpar@9: /* Column status */ alpar@9: #define DEAD_PRINCIPAL (-1) alpar@9: #define DEAD_NON_PRINCIPAL (-2) alpar@9: alpar@9: /* Macros for row and column status update and checking. */ alpar@9: #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) alpar@9: #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) alpar@9: #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) alpar@9: #define COL_IS_DEAD(c) (Col [c].start < ALIVE) alpar@9: #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) alpar@9: #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) alpar@9: #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } alpar@9: #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } alpar@9: #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Colamd reporting mechanism =========================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS) alpar@9: /* In MATLAB, matrices are 1-based to the user, but 0-based internally */ alpar@9: #define INDEX(i) ((i)+1) alpar@9: #else alpar@9: /* In C, matrices are 0-based and indices are reported as such in *_report */ alpar@9: #define INDEX(i) (i) alpar@9: #endif alpar@9: alpar@9: /* All output goes through the PRINTF macro. */ alpar@9: #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; } alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Prototypes of PRIVATE routines ======================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: PRIVATE Int init_rows_cols alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int p [], alpar@9: Int stats [COLAMD_STATS] alpar@9: ) ; alpar@9: alpar@9: PRIVATE void init_scoring alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int head [], alpar@9: double knobs [COLAMD_KNOBS], alpar@9: Int *p_n_row2, alpar@9: Int *p_n_col2, alpar@9: Int *p_max_deg alpar@9: ) ; alpar@9: alpar@9: PRIVATE Int find_ordering alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Int Alen, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int head [], alpar@9: Int n_col2, alpar@9: Int max_deg, alpar@9: Int pfree, alpar@9: Int aggressive alpar@9: ) ; alpar@9: alpar@9: PRIVATE void order_children alpar@9: ( alpar@9: Int n_col, alpar@9: Colamd_Col Col [], alpar@9: Int p [] alpar@9: ) ; alpar@9: alpar@9: PRIVATE void detect_super_cols alpar@9: ( alpar@9: alpar@9: #ifndef NDEBUG alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int head [], alpar@9: Int row_start, alpar@9: Int row_length alpar@9: ) ; alpar@9: alpar@9: PRIVATE Int garbage_collection alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int *pfree alpar@9: ) ; alpar@9: alpar@9: PRIVATE Int clear_mark alpar@9: ( alpar@9: Int tag_mark, alpar@9: Int max_mark, alpar@9: Int n_row, alpar@9: Colamd_Row Row [] alpar@9: ) ; alpar@9: alpar@9: PRIVATE void print_report alpar@9: ( alpar@9: char *method, alpar@9: Int stats [COLAMD_STATS] alpar@9: ) ; alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === Debugging prototypes and definitions ================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: alpar@9: #if 0 /* by mao */ alpar@9: #include alpar@9: #endif alpar@9: alpar@9: /* colamd_debug is the *ONLY* global variable, and is only */ alpar@9: /* present when debugging */ alpar@9: alpar@9: PRIVATE Int colamd_debug = 0 ; /* debug print level */ alpar@9: alpar@9: #define DEBUG0(params) { PRINTF (params) ; } alpar@9: #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; } alpar@9: #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; } alpar@9: #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; } alpar@9: #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; } alpar@9: alpar@9: #if 0 /* by mao */ alpar@9: #ifdef MATLAB_MEX_FILE alpar@9: #define ASSERT(expression) (mxAssert ((expression), "")) alpar@9: #else alpar@9: #define ASSERT(expression) (assert (expression)) alpar@9: #endif /* MATLAB_MEX_FILE */ alpar@9: #else alpar@9: #define ASSERT xassert alpar@9: #endif alpar@9: alpar@9: PRIVATE void colamd_get_debug /* gets the debug print level from getenv */ alpar@9: ( alpar@9: char *method alpar@9: ) ; alpar@9: alpar@9: PRIVATE void debug_deg_lists alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int head [], alpar@9: Int min_score, alpar@9: Int should, alpar@9: Int max_deg alpar@9: ) ; alpar@9: alpar@9: PRIVATE void debug_mark alpar@9: ( alpar@9: Int n_row, alpar@9: Colamd_Row Row [], alpar@9: Int tag_mark, alpar@9: Int max_mark alpar@9: ) ; alpar@9: alpar@9: PRIVATE void debug_matrix alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [] alpar@9: ) ; alpar@9: alpar@9: PRIVATE void debug_structures alpar@9: ( alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int n_col2 alpar@9: ) ; alpar@9: alpar@9: #else /* NDEBUG */ alpar@9: alpar@9: /* === No debugging ========================================================= */ alpar@9: alpar@9: #define DEBUG0(params) ; alpar@9: #define DEBUG1(params) ; alpar@9: #define DEBUG2(params) ; alpar@9: #define DEBUG3(params) ; alpar@9: #define DEBUG4(params) ; alpar@9: alpar@9: #define ASSERT(expression) alpar@9: alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === USER-CALLABLE ROUTINES: ============================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === colamd_recommended =================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: The colamd_recommended routine returns the suggested size for Alen. This alpar@9: value has been determined to provide good balance between the number of alpar@9: garbage collections and the memory requirements for colamd. If any alpar@9: argument is negative, or if integer overflow occurs, a 0 is returned as an alpar@9: error condition. 2*nnz space is required for the row and column alpar@9: indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is alpar@9: required for the Col and Row arrays, respectively, which are internal to alpar@9: colamd (roughly 6*n_col + 4*n_row). An additional n_col space is the alpar@9: minimal amount of "elbow room", and nnz/5 more space is recommended for alpar@9: run time efficiency. alpar@9: alpar@9: Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10. alpar@9: alpar@9: This function is not needed when using symamd. alpar@9: */ alpar@9: alpar@9: /* add two values of type size_t, and check for integer overflow */ alpar@9: static size_t t_add (size_t a, size_t b, int *ok) alpar@9: { alpar@9: (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ; alpar@9: return ((*ok) ? (a + b) : 0) ; alpar@9: } alpar@9: alpar@9: /* compute a*k where k is a small integer, and check for integer overflow */ alpar@9: static size_t t_mult (size_t a, size_t k, int *ok) alpar@9: { alpar@9: size_t i, s = 0 ; alpar@9: for (i = 0 ; i < k ; i++) alpar@9: { alpar@9: s = t_add (s, a, ok) ; alpar@9: } alpar@9: return (s) ; alpar@9: } alpar@9: alpar@9: /* size of the Col and Row structures */ alpar@9: #define COLAMD_C(n_col,ok) \ alpar@9: ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int))) alpar@9: alpar@9: #define COLAMD_R(n_row,ok) \ alpar@9: ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int))) alpar@9: alpar@9: alpar@9: PUBLIC size_t COLAMD_recommended /* returns recommended value of Alen. */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int nnz, /* number of nonzeros in A */ alpar@9: Int n_row, /* number of rows in A */ alpar@9: Int n_col /* number of columns in A */ alpar@9: ) alpar@9: { alpar@9: size_t s, c, r ; alpar@9: int ok = TRUE ; alpar@9: if (nnz < 0 || n_row < 0 || n_col < 0) alpar@9: { alpar@9: return (0) ; alpar@9: } alpar@9: s = t_mult (nnz, 2, &ok) ; /* 2*nnz */ alpar@9: c = COLAMD_C (n_col, &ok) ; /* size of column structures */ alpar@9: r = COLAMD_R (n_row, &ok) ; /* size of row structures */ alpar@9: s = t_add (s, c, &ok) ; alpar@9: s = t_add (s, r, &ok) ; alpar@9: s = t_add (s, n_col, &ok) ; /* elbow room */ alpar@9: s = t_add (s, nnz/5, &ok) ; /* elbow room */ alpar@9: ok = ok && (s < Int_MAX) ; alpar@9: return (ok ? s : 0) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === colamd_set_defaults ================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: The colamd_set_defaults routine sets the default values of the user- alpar@9: controllable parameters for colamd and symamd: alpar@9: alpar@9: Colamd: rows with more than max (16, knobs [0] * sqrt (n_col)) alpar@9: entries are removed prior to ordering. Columns with more than alpar@9: max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed alpar@9: prior to ordering, and placed last in the output column ordering. alpar@9: alpar@9: Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n)) alpar@9: entries are removed prior to ordering, and placed last in the alpar@9: output ordering. alpar@9: alpar@9: knobs [0] dense row control alpar@9: alpar@9: knobs [1] dense column control alpar@9: alpar@9: knobs [2] if nonzero, do aggresive absorption alpar@9: alpar@9: knobs [3..19] unused, but future versions might use this alpar@9: alpar@9: */ alpar@9: alpar@9: PUBLIC void COLAMD_set_defaults alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: double knobs [COLAMD_KNOBS] /* knob array */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int i ; alpar@9: alpar@9: if (!knobs) alpar@9: { alpar@9: return ; /* no knobs to initialize */ alpar@9: } alpar@9: for (i = 0 ; i < COLAMD_KNOBS ; i++) alpar@9: { alpar@9: knobs [i] = 0 ; alpar@9: } alpar@9: knobs [COLAMD_DENSE_ROW] = 10 ; alpar@9: knobs [COLAMD_DENSE_COL] = 10 ; alpar@9: knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/ alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === symamd =============================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: PUBLIC Int SYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n, /* number of rows and columns of A */ alpar@9: Int A [], /* row indices of A */ alpar@9: Int p [], /* column pointers of A */ alpar@9: Int perm [], /* output permutation, size n+1 */ alpar@9: double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */ alpar@9: Int stats [COLAMD_STATS], /* output statistics and error codes */ alpar@9: void * (*allocate) (size_t, size_t), alpar@9: /* pointer to calloc (ANSI C) or */ alpar@9: /* mxCalloc (for MATLAB mexFunction) */ alpar@9: void (*release) (void *) alpar@9: /* pointer to free (ANSI C) or */ alpar@9: /* mxFree (for MATLAB mexFunction) */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int *count ; /* length of each column of M, and col pointer*/ alpar@9: Int *mark ; /* mark array for finding duplicate entries */ alpar@9: Int *M ; /* row indices of matrix M */ alpar@9: size_t Mlen ; /* length of M */ alpar@9: Int n_row ; /* number of rows in M */ alpar@9: Int nnz ; /* number of entries in A */ alpar@9: Int i ; /* row index of A */ alpar@9: Int j ; /* column index of A */ alpar@9: Int k ; /* row index of M */ alpar@9: Int mnz ; /* number of nonzeros in M */ alpar@9: Int pp ; /* index into a column of A */ alpar@9: Int last_row ; /* last row seen in the current column */ alpar@9: Int length ; /* number of nonzeros in a column */ alpar@9: alpar@9: double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */ alpar@9: double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: colamd_get_debug ("symamd") ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Check the input arguments ======================================== */ alpar@9: alpar@9: if (!stats) alpar@9: { alpar@9: DEBUG0 (("symamd: stats not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: for (i = 0 ; i < COLAMD_STATS ; i++) alpar@9: { alpar@9: stats [i] = 0 ; alpar@9: } alpar@9: stats [COLAMD_STATUS] = COLAMD_OK ; alpar@9: stats [COLAMD_INFO1] = -1 ; alpar@9: stats [COLAMD_INFO2] = -1 ; alpar@9: alpar@9: if (!A) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; alpar@9: DEBUG0 (("symamd: A not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (!p) /* p is not present */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; alpar@9: DEBUG0 (("symamd: p not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (n < 0) /* n must be >= 0 */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; alpar@9: stats [COLAMD_INFO1] = n ; alpar@9: DEBUG0 (("symamd: n negative %d\n", n)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: nnz = p [n] ; alpar@9: if (nnz < 0) /* nnz must be >= 0 */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; alpar@9: stats [COLAMD_INFO1] = nnz ; alpar@9: DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (p [0] != 0) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; alpar@9: stats [COLAMD_INFO1] = p [0] ; alpar@9: DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: /* === If no knobs, set default knobs =================================== */ alpar@9: alpar@9: if (!knobs) alpar@9: { alpar@9: COLAMD_set_defaults (default_knobs) ; alpar@9: knobs = default_knobs ; alpar@9: } alpar@9: alpar@9: /* === Allocate count and mark ========================================== */ alpar@9: alpar@9: count = (Int *) ((*allocate) (n+1, sizeof (Int))) ; alpar@9: if (!count) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; alpar@9: DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ; alpar@9: if (!mark) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; alpar@9: (*release) ((void *) count) ; alpar@9: DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: /* === Compute column counts of M, check if A is valid ================== */ alpar@9: alpar@9: stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ alpar@9: alpar@9: for (i = 0 ; i < n ; i++) alpar@9: { alpar@9: mark [i] = -1 ; alpar@9: } alpar@9: alpar@9: for (j = 0 ; j < n ; j++) alpar@9: { alpar@9: last_row = -1 ; alpar@9: alpar@9: length = p [j+1] - p [j] ; alpar@9: if (length < 0) alpar@9: { alpar@9: /* column pointers must be non-decreasing */ alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; alpar@9: stats [COLAMD_INFO1] = j ; alpar@9: stats [COLAMD_INFO2] = length ; alpar@9: (*release) ((void *) count) ; alpar@9: (*release) ((void *) mark) ; alpar@9: DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: for (pp = p [j] ; pp < p [j+1] ; pp++) alpar@9: { alpar@9: i = A [pp] ; alpar@9: if (i < 0 || i >= n) alpar@9: { alpar@9: /* row index i, in column j, is out of bounds */ alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; alpar@9: stats [COLAMD_INFO1] = j ; alpar@9: stats [COLAMD_INFO2] = i ; alpar@9: stats [COLAMD_INFO3] = n ; alpar@9: (*release) ((void *) count) ; alpar@9: (*release) ((void *) mark) ; alpar@9: DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (i <= last_row || mark [i] == j) alpar@9: { alpar@9: /* row index is unsorted or repeated (or both), thus col */ alpar@9: /* is jumbled. This is a notice, not an error condition. */ alpar@9: stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; alpar@9: stats [COLAMD_INFO1] = j ; alpar@9: stats [COLAMD_INFO2] = i ; alpar@9: (stats [COLAMD_INFO3]) ++ ; alpar@9: DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ; alpar@9: } alpar@9: alpar@9: if (i > j && mark [i] != j) alpar@9: { alpar@9: /* row k of M will contain column indices i and j */ alpar@9: count [i]++ ; alpar@9: count [j]++ ; alpar@9: } alpar@9: alpar@9: /* mark the row as having been seen in this column */ alpar@9: mark [i] = j ; alpar@9: alpar@9: last_row = i ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* v2.4: removed free(mark) */ alpar@9: alpar@9: /* === Compute column pointers of M ===================================== */ alpar@9: alpar@9: /* use output permutation, perm, for column pointers of M */ alpar@9: perm [0] = 0 ; alpar@9: for (j = 1 ; j <= n ; j++) alpar@9: { alpar@9: perm [j] = perm [j-1] + count [j-1] ; alpar@9: } alpar@9: for (j = 0 ; j < n ; j++) alpar@9: { alpar@9: count [j] = perm [j] ; alpar@9: } alpar@9: alpar@9: /* === Construct M ====================================================== */ alpar@9: alpar@9: mnz = perm [n] ; alpar@9: n_row = mnz / 2 ; alpar@9: Mlen = COLAMD_recommended (mnz, n_row, n) ; alpar@9: M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ; alpar@9: DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n", alpar@9: n_row, n, mnz, (double) Mlen)) ; alpar@9: alpar@9: if (!M) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; alpar@9: (*release) ((void *) count) ; alpar@9: (*release) ((void *) mark) ; alpar@9: DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: k = 0 ; alpar@9: alpar@9: if (stats [COLAMD_STATUS] == COLAMD_OK) alpar@9: { alpar@9: /* Matrix is OK */ alpar@9: for (j = 0 ; j < n ; j++) alpar@9: { alpar@9: ASSERT (p [j+1] - p [j] >= 0) ; alpar@9: for (pp = p [j] ; pp < p [j+1] ; pp++) alpar@9: { alpar@9: i = A [pp] ; alpar@9: ASSERT (i >= 0 && i < n) ; alpar@9: if (i > j) alpar@9: { alpar@9: /* row k of M contains column indices i and j */ alpar@9: M [count [i]++] = k ; alpar@9: M [count [j]++] = k ; alpar@9: k++ ; alpar@9: } alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: { alpar@9: /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */ alpar@9: DEBUG0 (("symamd: Duplicates in A.\n")) ; alpar@9: for (i = 0 ; i < n ; i++) alpar@9: { alpar@9: mark [i] = -1 ; alpar@9: } alpar@9: for (j = 0 ; j < n ; j++) alpar@9: { alpar@9: ASSERT (p [j+1] - p [j] >= 0) ; alpar@9: for (pp = p [j] ; pp < p [j+1] ; pp++) alpar@9: { alpar@9: i = A [pp] ; alpar@9: ASSERT (i >= 0 && i < n) ; alpar@9: if (i > j && mark [i] != j) alpar@9: { alpar@9: /* row k of M contains column indices i and j */ alpar@9: M [count [i]++] = k ; alpar@9: M [count [j]++] = k ; alpar@9: k++ ; alpar@9: mark [i] = j ; alpar@9: } alpar@9: } alpar@9: } alpar@9: /* v2.4: free(mark) moved below */ alpar@9: } alpar@9: alpar@9: /* count and mark no longer needed */ alpar@9: (*release) ((void *) count) ; alpar@9: (*release) ((void *) mark) ; /* v2.4: free (mark) moved here */ alpar@9: ASSERT (k == n_row) ; alpar@9: alpar@9: /* === Adjust the knobs for M =========================================== */ alpar@9: alpar@9: for (i = 0 ; i < COLAMD_KNOBS ; i++) alpar@9: { alpar@9: cknobs [i] = knobs [i] ; alpar@9: } alpar@9: alpar@9: /* there are no dense rows in M */ alpar@9: cknobs [COLAMD_DENSE_ROW] = -1 ; alpar@9: cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ; alpar@9: alpar@9: /* === Order the columns of M =========================================== */ alpar@9: alpar@9: /* v2.4: colamd cannot fail here, so the error check is removed */ alpar@9: (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ; alpar@9: alpar@9: /* Note that the output permutation is now in perm */ alpar@9: alpar@9: /* === get the statistics for symamd from colamd ======================== */ alpar@9: alpar@9: /* a dense column in colamd means a dense row and col in symamd */ alpar@9: stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ; alpar@9: alpar@9: /* === Free M =========================================================== */ alpar@9: alpar@9: (*release) ((void *) M) ; alpar@9: DEBUG0 (("symamd: done.\n")) ; alpar@9: return (TRUE) ; alpar@9: alpar@9: } alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === colamd =============================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: The colamd routine computes a column ordering Q of a sparse matrix alpar@9: A such that the LU factorization P(AQ) = LU remains sparse, where P is alpar@9: selected via partial pivoting. The routine can also be viewed as alpar@9: providing a permutation Q such that the Cholesky factorization alpar@9: (AQ)'(AQ) = LL' remains sparse. alpar@9: */ alpar@9: alpar@9: PUBLIC Int COLAMD_MAIN /* returns TRUE if successful, FALSE otherwise*/ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, /* number of rows in A */ alpar@9: Int n_col, /* number of columns in A */ alpar@9: Int Alen, /* length of A */ alpar@9: Int A [], /* row indices of A */ alpar@9: Int p [], /* pointers to columns in A */ alpar@9: double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */ alpar@9: Int stats [COLAMD_STATS] /* output statistics and error codes */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int i ; /* loop index */ alpar@9: Int nnz ; /* nonzeros in A */ alpar@9: size_t Row_size ; /* size of Row [], in integers */ alpar@9: size_t Col_size ; /* size of Col [], in integers */ alpar@9: size_t need ; /* minimum required length of A */ alpar@9: Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ alpar@9: Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ alpar@9: Int n_col2 ; /* number of non-dense, non-empty columns */ alpar@9: Int n_row2 ; /* number of non-dense, non-empty rows */ alpar@9: Int ngarbage ; /* number of garbage collections performed */ alpar@9: Int max_deg ; /* maximum row degree */ alpar@9: double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ alpar@9: Int aggressive ; /* do aggressive absorption */ alpar@9: int ok ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: colamd_get_debug ("colamd") ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Check the input arguments ======================================== */ alpar@9: alpar@9: if (!stats) alpar@9: { alpar@9: DEBUG0 (("colamd: stats not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: for (i = 0 ; i < COLAMD_STATS ; i++) alpar@9: { alpar@9: stats [i] = 0 ; alpar@9: } alpar@9: stats [COLAMD_STATUS] = COLAMD_OK ; alpar@9: stats [COLAMD_INFO1] = -1 ; alpar@9: stats [COLAMD_INFO2] = -1 ; alpar@9: alpar@9: if (!A) /* A is not present */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; alpar@9: DEBUG0 (("colamd: A not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (!p) /* p is not present */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; alpar@9: DEBUG0 (("colamd: p not present\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (n_row < 0) /* n_row must be >= 0 */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; alpar@9: stats [COLAMD_INFO1] = n_row ; alpar@9: DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (n_col < 0) /* n_col must be >= 0 */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; alpar@9: stats [COLAMD_INFO1] = n_col ; alpar@9: DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: nnz = p [n_col] ; alpar@9: if (nnz < 0) /* nnz must be >= 0 */ alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; alpar@9: stats [COLAMD_INFO1] = nnz ; alpar@9: DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (p [0] != 0) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; alpar@9: stats [COLAMD_INFO1] = p [0] ; alpar@9: DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: /* === If no knobs, set default knobs =================================== */ alpar@9: alpar@9: if (!knobs) alpar@9: { alpar@9: COLAMD_set_defaults (default_knobs) ; alpar@9: knobs = default_knobs ; alpar@9: } alpar@9: alpar@9: aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ; alpar@9: alpar@9: /* === Allocate the Row and Col arrays from array A ===================== */ alpar@9: alpar@9: ok = TRUE ; alpar@9: Col_size = COLAMD_C (n_col, &ok) ; /* size of Col array of structs */ alpar@9: Row_size = COLAMD_R (n_row, &ok) ; /* size of Row array of structs */ alpar@9: alpar@9: /* need = 2*nnz + n_col + Col_size + Row_size ; */ alpar@9: need = t_mult (nnz, 2, &ok) ; alpar@9: need = t_add (need, n_col, &ok) ; alpar@9: need = t_add (need, Col_size, &ok) ; alpar@9: need = t_add (need, Row_size, &ok) ; alpar@9: alpar@9: if (!ok || need > (size_t) Alen || need > Int_MAX) alpar@9: { alpar@9: /* not enough space in array A to perform the ordering */ alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; alpar@9: stats [COLAMD_INFO1] = need ; alpar@9: stats [COLAMD_INFO2] = Alen ; alpar@9: DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: Alen -= Col_size + Row_size ; alpar@9: Col = (Colamd_Col *) &A [Alen] ; alpar@9: Row = (Colamd_Row *) &A [Alen + Col_size] ; alpar@9: alpar@9: /* === Construct the row and column data structures ===================== */ alpar@9: alpar@9: if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) alpar@9: { alpar@9: /* input matrix is invalid */ alpar@9: DEBUG0 (("colamd: Matrix invalid\n")) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: /* === Initialize scores, kill dense rows/columns ======================= */ alpar@9: alpar@9: init_scoring (n_row, n_col, Row, Col, A, p, knobs, alpar@9: &n_row2, &n_col2, &max_deg) ; alpar@9: alpar@9: /* === Order the supercolumns =========================================== */ alpar@9: alpar@9: ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, alpar@9: n_col2, max_deg, 2*nnz, aggressive) ; alpar@9: alpar@9: /* === Order the non-principal columns ================================== */ alpar@9: alpar@9: order_children (n_col, Col, p) ; alpar@9: alpar@9: /* === Return statistics in stats ======================================= */ alpar@9: alpar@9: stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; alpar@9: stats [COLAMD_DENSE_COL] = n_col - n_col2 ; alpar@9: stats [COLAMD_DEFRAG_COUNT] = ngarbage ; alpar@9: DEBUG0 (("colamd: done.\n")) ; alpar@9: return (TRUE) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === colamd_report ======================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: PUBLIC void COLAMD_report alpar@9: ( alpar@9: Int stats [COLAMD_STATS] alpar@9: ) alpar@9: { alpar@9: print_report ("colamd", stats) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === symamd_report ======================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: PUBLIC void SYMAMD_report alpar@9: ( alpar@9: Int stats [COLAMD_STATS] alpar@9: ) alpar@9: { alpar@9: print_report ("symamd", stats) ; alpar@9: } alpar@9: alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === NON-USER-CALLABLE ROUTINES: ========================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* There are no user-callable routines beyond this point in the file */ alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === init_rows_cols ======================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Takes the column form of the matrix in A and creates the row form of the alpar@9: matrix. Also, row and column attributes are stored in the Col and Row alpar@9: structs. If the columns are un-sorted or contain duplicate row indices, alpar@9: this routine will also sort and remove duplicate row indices from the alpar@9: column form of the matrix. Returns FALSE if the matrix is invalid, alpar@9: TRUE otherwise. Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, /* number of rows of A */ alpar@9: Int n_col, /* number of columns of A */ alpar@9: Colamd_Row Row [], /* of size n_row+1 */ alpar@9: Colamd_Col Col [], /* of size n_col+1 */ alpar@9: Int A [], /* row indices of A, of size Alen */ alpar@9: Int p [], /* pointers to columns in A, of size n_col+1 */ alpar@9: Int stats [COLAMD_STATS] /* colamd statistics */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int col ; /* a column index */ alpar@9: Int row ; /* a row index */ alpar@9: Int *cp ; /* a column pointer */ alpar@9: Int *cp_end ; /* a pointer to the end of a column */ alpar@9: Int *rp ; /* a row pointer */ alpar@9: Int *rp_end ; /* a pointer to the end of a row */ alpar@9: Int last_row ; /* previous row */ alpar@9: alpar@9: /* === Initialize columns, and check column pointers ==================== */ alpar@9: alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: Col [col].start = p [col] ; alpar@9: Col [col].length = p [col+1] - p [col] ; alpar@9: alpar@9: if (Col [col].length < 0) alpar@9: { alpar@9: /* column pointers must be non-decreasing */ alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; alpar@9: stats [COLAMD_INFO1] = col ; alpar@9: stats [COLAMD_INFO2] = Col [col].length ; alpar@9: DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: Col [col].shared1.thickness = 1 ; alpar@9: Col [col].shared2.score = 0 ; alpar@9: Col [col].shared3.prev = EMPTY ; alpar@9: Col [col].shared4.degree_next = EMPTY ; alpar@9: } alpar@9: alpar@9: /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ alpar@9: alpar@9: /* === Scan columns, compute row degrees, and check row indices ========= */ alpar@9: alpar@9: stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ alpar@9: alpar@9: for (row = 0 ; row < n_row ; row++) alpar@9: { alpar@9: Row [row].length = 0 ; alpar@9: Row [row].shared2.mark = -1 ; alpar@9: } alpar@9: alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: last_row = -1 ; alpar@9: alpar@9: cp = &A [p [col]] ; alpar@9: cp_end = &A [p [col+1]] ; alpar@9: alpar@9: while (cp < cp_end) alpar@9: { alpar@9: row = *cp++ ; alpar@9: alpar@9: /* make sure row indices within range */ alpar@9: if (row < 0 || row >= n_row) alpar@9: { alpar@9: stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; alpar@9: stats [COLAMD_INFO1] = col ; alpar@9: stats [COLAMD_INFO2] = row ; alpar@9: stats [COLAMD_INFO3] = n_row ; alpar@9: DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; alpar@9: return (FALSE) ; alpar@9: } alpar@9: alpar@9: if (row <= last_row || Row [row].shared2.mark == col) alpar@9: { alpar@9: /* row index are unsorted or repeated (or both), thus col */ alpar@9: /* is jumbled. This is a notice, not an error condition. */ alpar@9: stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; alpar@9: stats [COLAMD_INFO1] = col ; alpar@9: stats [COLAMD_INFO2] = row ; alpar@9: (stats [COLAMD_INFO3]) ++ ; alpar@9: DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); alpar@9: } alpar@9: alpar@9: if (Row [row].shared2.mark != col) alpar@9: { alpar@9: Row [row].length++ ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* this is a repeated entry in the column, */ alpar@9: /* it will be removed */ alpar@9: Col [col].length-- ; alpar@9: } alpar@9: alpar@9: /* mark the row as having been seen in this column */ alpar@9: Row [row].shared2.mark = col ; alpar@9: alpar@9: last_row = row ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Compute row pointers ============================================= */ alpar@9: alpar@9: /* row form of the matrix starts directly after the column */ alpar@9: /* form of matrix in A */ alpar@9: Row [0].start = p [n_col] ; alpar@9: Row [0].shared1.p = Row [0].start ; alpar@9: Row [0].shared2.mark = -1 ; alpar@9: for (row = 1 ; row < n_row ; row++) alpar@9: { alpar@9: Row [row].start = Row [row-1].start + Row [row-1].length ; alpar@9: Row [row].shared1.p = Row [row].start ; alpar@9: Row [row].shared2.mark = -1 ; alpar@9: } alpar@9: alpar@9: /* === Create row form ================================================== */ alpar@9: alpar@9: if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) alpar@9: { alpar@9: /* if cols jumbled, watch for repeated row indices */ alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: cp = &A [p [col]] ; alpar@9: cp_end = &A [p [col+1]] ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: row = *cp++ ; alpar@9: if (Row [row].shared2.mark != col) alpar@9: { alpar@9: A [(Row [row].shared1.p)++] = col ; alpar@9: Row [row].shared2.mark = col ; alpar@9: } alpar@9: } alpar@9: } alpar@9: } alpar@9: else alpar@9: { alpar@9: /* if cols not jumbled, we don't need the mark (this is faster) */ alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: cp = &A [p [col]] ; alpar@9: cp_end = &A [p [col+1]] ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: A [(Row [*cp++].shared1.p)++] = col ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Clear the row marks and set row degrees ========================== */ alpar@9: alpar@9: for (row = 0 ; row < n_row ; row++) alpar@9: { alpar@9: Row [row].shared2.mark = 0 ; alpar@9: Row [row].shared1.degree = Row [row].length ; alpar@9: } alpar@9: alpar@9: /* === See if we need to re-create columns ============================== */ alpar@9: alpar@9: if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) alpar@9: { alpar@9: DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: /* make sure column lengths are correct */ alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: p [col] = Col [col].length ; alpar@9: } alpar@9: for (row = 0 ; row < n_row ; row++) alpar@9: { alpar@9: rp = &A [Row [row].start] ; alpar@9: rp_end = rp + Row [row].length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: p [*rp++]-- ; alpar@9: } alpar@9: } alpar@9: for (col = 0 ; col < n_col ; col++) alpar@9: { alpar@9: ASSERT (p [col] == 0) ; alpar@9: } alpar@9: /* now p is all zero (different than when debugging is turned off) */ alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Compute col pointers ========================================= */ alpar@9: alpar@9: /* col form of the matrix starts at A [0]. */ alpar@9: /* Note, we may have a gap between the col form and the row */ alpar@9: /* form if there were duplicate entries, if so, it will be */ alpar@9: /* removed upon the first garbage collection */ alpar@9: Col [0].start = 0 ; alpar@9: p [0] = Col [0].start ; alpar@9: for (col = 1 ; col < n_col ; col++) alpar@9: { alpar@9: /* note that the lengths here are for pruned columns, i.e. */ alpar@9: /* no duplicate row indices will exist for these columns */ alpar@9: Col [col].start = Col [col-1].start + Col [col-1].length ; alpar@9: p [col] = Col [col].start ; alpar@9: } alpar@9: alpar@9: /* === Re-create col form =========================================== */ alpar@9: alpar@9: for (row = 0 ; row < n_row ; row++) alpar@9: { alpar@9: rp = &A [Row [row].start] ; alpar@9: rp_end = rp + Row [row].length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: A [(p [*rp++])++] = row ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Done. Matrix is not (or no longer) jumbled ====================== */ alpar@9: alpar@9: return (TRUE) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === init_scoring ========================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Kills dense or empty columns and rows, calculates an initial score for alpar@9: each column, and places all columns in the degree lists. Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE void init_scoring alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, /* number of rows of A */ alpar@9: Int n_col, /* number of columns of A */ alpar@9: Colamd_Row Row [], /* of size n_row+1 */ alpar@9: Colamd_Col Col [], /* of size n_col+1 */ alpar@9: Int A [], /* column form and row form of A */ alpar@9: Int head [], /* of size n_col+1 */ alpar@9: double knobs [COLAMD_KNOBS],/* parameters */ alpar@9: Int *p_n_row2, /* number of non-dense, non-empty rows */ alpar@9: Int *p_n_col2, /* number of non-dense, non-empty columns */ alpar@9: Int *p_max_deg /* maximum row degree */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int c ; /* a column index */ alpar@9: Int r, row ; /* a row index */ alpar@9: Int *cp ; /* a column pointer */ alpar@9: Int deg ; /* degree of a row or column */ alpar@9: Int *cp_end ; /* a pointer to the end of a column */ alpar@9: Int *new_cp ; /* new column pointer */ alpar@9: Int col_length ; /* length of pruned column */ alpar@9: Int score ; /* current column score */ alpar@9: Int n_col2 ; /* number of non-dense, non-empty columns */ alpar@9: Int n_row2 ; /* number of non-dense, non-empty rows */ alpar@9: Int dense_row_count ; /* remove rows with more entries than this */ alpar@9: Int dense_col_count ; /* remove cols with more entries than this */ alpar@9: Int min_score ; /* smallest column score */ alpar@9: Int max_deg ; /* maximum row degree */ alpar@9: Int next_col ; /* Used to add to degree list.*/ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: Int debug_count ; /* debug only. */ alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Extract knobs ==================================================== */ alpar@9: alpar@9: /* Note: if knobs contains a NaN, this is undefined: */ alpar@9: if (knobs [COLAMD_DENSE_ROW] < 0) alpar@9: { alpar@9: /* only remove completely dense rows */ alpar@9: dense_row_count = n_col-1 ; alpar@9: } alpar@9: else alpar@9: { alpar@9: dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ; alpar@9: } alpar@9: if (knobs [COLAMD_DENSE_COL] < 0) alpar@9: { alpar@9: /* only remove completely dense columns */ alpar@9: dense_col_count = n_row-1 ; alpar@9: } alpar@9: else alpar@9: { alpar@9: dense_col_count = alpar@9: DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ; alpar@9: } alpar@9: alpar@9: DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; alpar@9: max_deg = 0 ; alpar@9: n_col2 = n_col ; alpar@9: n_row2 = n_row ; alpar@9: alpar@9: /* === Kill empty columns =============================================== */ alpar@9: alpar@9: /* Put the empty columns at the end in their natural order, so that LU */ alpar@9: /* factorization can proceed as far as possible. */ alpar@9: for (c = n_col-1 ; c >= 0 ; c--) alpar@9: { alpar@9: deg = Col [c].length ; alpar@9: if (deg == 0) alpar@9: { alpar@9: /* this is a empty column, kill and order it last */ alpar@9: Col [c].shared2.order = --n_col2 ; alpar@9: KILL_PRINCIPAL_COL (c) ; alpar@9: } alpar@9: } alpar@9: DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; alpar@9: alpar@9: /* === Kill dense columns =============================================== */ alpar@9: alpar@9: /* Put the dense columns at the end, in their natural order */ alpar@9: for (c = n_col-1 ; c >= 0 ; c--) alpar@9: { alpar@9: /* skip any dead columns */ alpar@9: if (COL_IS_DEAD (c)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: deg = Col [c].length ; alpar@9: if (deg > dense_col_count) alpar@9: { alpar@9: /* this is a dense column, kill and order it last */ alpar@9: Col [c].shared2.order = --n_col2 ; alpar@9: /* decrement the row degrees */ alpar@9: cp = &A [Col [c].start] ; alpar@9: cp_end = cp + Col [c].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: Row [*cp++].shared1.degree-- ; alpar@9: } alpar@9: KILL_PRINCIPAL_COL (c) ; alpar@9: } alpar@9: } alpar@9: DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; alpar@9: alpar@9: /* === Kill dense and empty rows ======================================== */ alpar@9: alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: deg = Row [r].shared1.degree ; alpar@9: ASSERT (deg >= 0 && deg <= n_col) ; alpar@9: if (deg > dense_row_count || deg == 0) alpar@9: { alpar@9: /* kill a dense or empty row */ alpar@9: KILL_ROW (r) ; alpar@9: --n_row2 ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* keep track of max degree of remaining rows */ alpar@9: max_deg = MAX (max_deg, deg) ; alpar@9: } alpar@9: } alpar@9: DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; alpar@9: alpar@9: /* === Compute initial column scores ==================================== */ alpar@9: alpar@9: /* At this point the row degrees are accurate. They reflect the number */ alpar@9: /* of "live" (non-dense) columns in each row. No empty rows exist. */ alpar@9: /* Some "live" columns may contain only dead rows, however. These are */ alpar@9: /* pruned in the code below. */ alpar@9: alpar@9: /* now find the initial matlab score for each column */ alpar@9: for (c = n_col-1 ; c >= 0 ; c--) alpar@9: { alpar@9: /* skip dead column */ alpar@9: if (COL_IS_DEAD (c)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: score = 0 ; alpar@9: cp = &A [Col [c].start] ; alpar@9: new_cp = cp ; alpar@9: cp_end = cp + Col [c].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: /* get a row */ alpar@9: row = *cp++ ; alpar@9: /* skip if dead */ alpar@9: if (ROW_IS_DEAD (row)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: /* compact the column */ alpar@9: *new_cp++ = row ; alpar@9: /* add row's external degree */ alpar@9: score += Row [row].shared1.degree - 1 ; alpar@9: /* guard against integer overflow */ alpar@9: score = MIN (score, n_col) ; alpar@9: } alpar@9: /* determine pruned column length */ alpar@9: col_length = (Int) (new_cp - &A [Col [c].start]) ; alpar@9: if (col_length == 0) alpar@9: { alpar@9: /* a newly-made null column (all rows in this col are "dense" */ alpar@9: /* and have already been killed) */ alpar@9: DEBUG2 (("Newly null killed: %d\n", c)) ; alpar@9: Col [c].shared2.order = --n_col2 ; alpar@9: KILL_PRINCIPAL_COL (c) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* set column length and set score */ alpar@9: ASSERT (score >= 0) ; alpar@9: ASSERT (score <= n_col) ; alpar@9: Col [c].length = col_length ; alpar@9: Col [c].shared2.score = score ; alpar@9: } alpar@9: } alpar@9: DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", alpar@9: n_col-n_col2)) ; alpar@9: alpar@9: /* At this point, all empty rows and columns are dead. All live columns */ alpar@9: /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ alpar@9: /* yet). Rows may contain dead columns, but all live rows contain at */ alpar@9: /* least one live column. */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_structures (n_row, n_col, Row, Col, A, n_col2) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Initialize degree lists ========================================== */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_count = 0 ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* clear the hash buckets */ alpar@9: for (c = 0 ; c <= n_col ; c++) alpar@9: { alpar@9: head [c] = EMPTY ; alpar@9: } alpar@9: min_score = n_col ; alpar@9: /* place in reverse order, so low column indices are at the front */ alpar@9: /* of the lists. This is to encourage natural tie-breaking */ alpar@9: for (c = n_col-1 ; c >= 0 ; c--) alpar@9: { alpar@9: /* only add principal columns to degree lists */ alpar@9: if (COL_IS_ALIVE (c)) alpar@9: { alpar@9: DEBUG4 (("place %d score %d minscore %d ncol %d\n", alpar@9: c, Col [c].shared2.score, min_score, n_col)) ; alpar@9: alpar@9: /* === Add columns score to DList =============================== */ alpar@9: alpar@9: score = Col [c].shared2.score ; alpar@9: alpar@9: ASSERT (min_score >= 0) ; alpar@9: ASSERT (min_score <= n_col) ; alpar@9: ASSERT (score >= 0) ; alpar@9: ASSERT (score <= n_col) ; alpar@9: ASSERT (head [score] >= EMPTY) ; alpar@9: alpar@9: /* now add this column to dList at proper score location */ alpar@9: next_col = head [score] ; alpar@9: Col [c].shared3.prev = EMPTY ; alpar@9: Col [c].shared4.degree_next = next_col ; alpar@9: alpar@9: /* if there already was a column with the same score, set its */ alpar@9: /* previous pointer to this new column */ alpar@9: if (next_col != EMPTY) alpar@9: { alpar@9: Col [next_col].shared3.prev = c ; alpar@9: } alpar@9: head [score] = c ; alpar@9: alpar@9: /* see if this score is less than current min */ alpar@9: min_score = MIN (min_score, score) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_count++ ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: } alpar@9: } alpar@9: alpar@9: #ifndef NDEBUG alpar@9: DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n", alpar@9: debug_count, n_col, n_col-debug_count)) ; alpar@9: ASSERT (debug_count == n_col2) ; alpar@9: debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Return number of remaining columns, and max row degree =========== */ alpar@9: alpar@9: *p_n_col2 = n_col2 ; alpar@9: *p_n_row2 = n_row2 ; alpar@9: *p_max_deg = max_deg ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === find_ordering ======================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Order the principal columns of the supercolumn form of the matrix alpar@9: (no supercolumns on input). Uses a minimum approximate column minimum alpar@9: degree ordering method. Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE Int find_ordering /* return the number of garbage collections */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, /* number of rows of A */ alpar@9: Int n_col, /* number of columns of A */ alpar@9: Int Alen, /* size of A, 2*nnz + n_col or larger */ alpar@9: Colamd_Row Row [], /* of size n_row+1 */ alpar@9: Colamd_Col Col [], /* of size n_col+1 */ alpar@9: Int A [], /* column form and row form of A */ alpar@9: Int head [], /* of size n_col+1 */ alpar@9: Int n_col2, /* Remaining columns to order */ alpar@9: Int max_deg, /* Maximum row degree */ alpar@9: Int pfree, /* index of first free slot (2*nnz on entry) */ alpar@9: Int aggressive alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int k ; /* current pivot ordering step */ alpar@9: Int pivot_col ; /* current pivot column */ alpar@9: Int *cp ; /* a column pointer */ alpar@9: Int *rp ; /* a row pointer */ alpar@9: Int pivot_row ; /* current pivot row */ alpar@9: Int *new_cp ; /* modified column pointer */ alpar@9: Int *new_rp ; /* modified row pointer */ alpar@9: Int pivot_row_start ; /* pointer to start of pivot row */ alpar@9: Int pivot_row_degree ; /* number of columns in pivot row */ alpar@9: Int pivot_row_length ; /* number of supercolumns in pivot row */ alpar@9: Int pivot_col_score ; /* score of pivot column */ alpar@9: Int needed_memory ; /* free space needed for pivot row */ alpar@9: Int *cp_end ; /* pointer to the end of a column */ alpar@9: Int *rp_end ; /* pointer to the end of a row */ alpar@9: Int row ; /* a row index */ alpar@9: Int col ; /* a column index */ alpar@9: Int max_score ; /* maximum possible score */ alpar@9: Int cur_score ; /* score of current column */ alpar@9: unsigned Int hash ; /* hash value for supernode detection */ alpar@9: Int head_column ; /* head of hash bucket */ alpar@9: Int first_col ; /* first column in hash bucket */ alpar@9: Int tag_mark ; /* marker value for mark array */ alpar@9: Int row_mark ; /* Row [row].shared2.mark */ alpar@9: Int set_difference ; /* set difference size of row with pivot row */ alpar@9: Int min_score ; /* smallest column score */ alpar@9: Int col_thickness ; /* "thickness" (no. of columns in a supercol) */ alpar@9: Int max_mark ; /* maximum value of tag_mark */ alpar@9: Int pivot_col_thickness ; /* number of columns represented by pivot col */ alpar@9: Int prev_col ; /* Used by Dlist operations. */ alpar@9: Int next_col ; /* Used by Dlist operations. */ alpar@9: Int ngarbage ; /* number of garbage collections performed */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: Int debug_d ; /* debug loop counter */ alpar@9: Int debug_step = 0 ; /* debug loop counter */ alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Initialization and clear mark ==================================== */ alpar@9: alpar@9: max_mark = INT_MAX - n_col ; /* INT_MAX defined in */ alpar@9: tag_mark = clear_mark (0, max_mark, n_row, Row) ; alpar@9: min_score = 0 ; alpar@9: ngarbage = 0 ; alpar@9: DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; alpar@9: alpar@9: /* === Order the columns ================================================ */ alpar@9: alpar@9: for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) alpar@9: { alpar@9: alpar@9: #ifndef NDEBUG alpar@9: if (debug_step % 100 == 0) alpar@9: { alpar@9: DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; alpar@9: } alpar@9: debug_step++ ; alpar@9: debug_deg_lists (n_row, n_col, Row, Col, head, alpar@9: min_score, n_col2-k, max_deg) ; alpar@9: debug_matrix (n_row, n_col, Row, Col, A) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Select pivot column, and order it ============================ */ alpar@9: alpar@9: /* make sure degree list isn't empty */ alpar@9: ASSERT (min_score >= 0) ; alpar@9: ASSERT (min_score <= n_col) ; alpar@9: ASSERT (head [min_score] >= EMPTY) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: for (debug_d = 0 ; debug_d < min_score ; debug_d++) alpar@9: { alpar@9: ASSERT (head [debug_d] == EMPTY) ; alpar@9: } alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* get pivot column from head of minimum degree list */ alpar@9: while (head [min_score] == EMPTY && min_score < n_col) alpar@9: { alpar@9: min_score++ ; alpar@9: } alpar@9: pivot_col = head [min_score] ; alpar@9: ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; alpar@9: next_col = Col [pivot_col].shared4.degree_next ; alpar@9: head [min_score] = next_col ; alpar@9: if (next_col != EMPTY) alpar@9: { alpar@9: Col [next_col].shared3.prev = EMPTY ; alpar@9: } alpar@9: alpar@9: ASSERT (COL_IS_ALIVE (pivot_col)) ; alpar@9: alpar@9: /* remember score for defrag check */ alpar@9: pivot_col_score = Col [pivot_col].shared2.score ; alpar@9: alpar@9: /* the pivot column is the kth column in the pivot order */ alpar@9: Col [pivot_col].shared2.order = k ; alpar@9: alpar@9: /* increment order count by column thickness */ alpar@9: pivot_col_thickness = Col [pivot_col].shared1.thickness ; alpar@9: k += pivot_col_thickness ; alpar@9: ASSERT (pivot_col_thickness > 0) ; alpar@9: DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ; alpar@9: alpar@9: /* === Garbage_collection, if necessary ============================= */ alpar@9: alpar@9: needed_memory = MIN (pivot_col_score, n_col - k) ; alpar@9: if (pfree + needed_memory >= Alen) alpar@9: { alpar@9: pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; alpar@9: ngarbage++ ; alpar@9: /* after garbage collection we will have enough */ alpar@9: ASSERT (pfree + needed_memory < Alen) ; alpar@9: /* garbage collection has wiped out the Row[].shared2.mark array */ alpar@9: tag_mark = clear_mark (0, max_mark, n_row, Row) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_matrix (n_row, n_col, Row, Col, A) ; alpar@9: #endif /* NDEBUG */ alpar@9: } alpar@9: alpar@9: /* === Compute pivot row pattern ==================================== */ alpar@9: alpar@9: /* get starting location for this new merged row */ alpar@9: pivot_row_start = pfree ; alpar@9: alpar@9: /* initialize new row counts to zero */ alpar@9: pivot_row_degree = 0 ; alpar@9: alpar@9: /* tag pivot column as having been visited so it isn't included */ alpar@9: /* in merged pivot row */ alpar@9: Col [pivot_col].shared1.thickness = -pivot_col_thickness ; alpar@9: alpar@9: /* pivot row is the union of all rows in the pivot column pattern */ alpar@9: cp = &A [Col [pivot_col].start] ; alpar@9: cp_end = cp + Col [pivot_col].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: /* get a row */ alpar@9: row = *cp++ ; alpar@9: DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; alpar@9: /* skip if row is dead */ alpar@9: if (ROW_IS_ALIVE (row)) alpar@9: { alpar@9: rp = &A [Row [row].start] ; alpar@9: rp_end = rp + Row [row].length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: /* get a column */ alpar@9: col = *rp++ ; alpar@9: /* add the column, if alive and untagged */ alpar@9: col_thickness = Col [col].shared1.thickness ; alpar@9: if (col_thickness > 0 && COL_IS_ALIVE (col)) alpar@9: { alpar@9: /* tag column in pivot row */ alpar@9: Col [col].shared1.thickness = -col_thickness ; alpar@9: ASSERT (pfree < Alen) ; alpar@9: /* place column in pivot row */ alpar@9: A [pfree++] = col ; alpar@9: pivot_row_degree += col_thickness ; alpar@9: } alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: /* clear tag on pivot column */ alpar@9: Col [pivot_col].shared1.thickness = pivot_col_thickness ; alpar@9: max_deg = MAX (max_deg, pivot_row_degree) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: DEBUG3 (("check2\n")) ; alpar@9: debug_mark (n_row, Row, tag_mark, max_mark) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Kill all rows used to construct pivot row ==================== */ alpar@9: alpar@9: /* also kill pivot row, temporarily */ alpar@9: cp = &A [Col [pivot_col].start] ; alpar@9: cp_end = cp + Col [pivot_col].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: /* may be killing an already dead row */ alpar@9: row = *cp++ ; alpar@9: DEBUG3 (("Kill row in pivot col: %d\n", row)) ; alpar@9: KILL_ROW (row) ; alpar@9: } alpar@9: alpar@9: /* === Select a row index to use as the new pivot row =============== */ alpar@9: alpar@9: pivot_row_length = pfree - pivot_row_start ; alpar@9: if (pivot_row_length > 0) alpar@9: { alpar@9: /* pick the "pivot" row arbitrarily (first row in col) */ alpar@9: pivot_row = A [Col [pivot_col].start] ; alpar@9: DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* there is no pivot row, since it is of zero length */ alpar@9: pivot_row = EMPTY ; alpar@9: ASSERT (pivot_row_length == 0) ; alpar@9: } alpar@9: ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; alpar@9: alpar@9: /* === Approximate degree computation =============================== */ alpar@9: alpar@9: /* Here begins the computation of the approximate degree. The column */ alpar@9: /* score is the sum of the pivot row "length", plus the size of the */ alpar@9: /* set differences of each row in the column minus the pattern of the */ alpar@9: /* pivot row itself. The column ("thickness") itself is also */ alpar@9: /* excluded from the column score (we thus use an approximate */ alpar@9: /* external degree). */ alpar@9: alpar@9: /* The time taken by the following code (compute set differences, and */ alpar@9: /* add them up) is proportional to the size of the data structure */ alpar@9: /* being scanned - that is, the sum of the sizes of each column in */ alpar@9: /* the pivot row. Thus, the amortized time to compute a column score */ alpar@9: /* is proportional to the size of that column (where size, in this */ alpar@9: /* context, is the column "length", or the number of row indices */ alpar@9: /* in that column). The number of row indices in a column is */ alpar@9: /* monotonically non-decreasing, from the length of the original */ alpar@9: /* column on input to colamd. */ alpar@9: alpar@9: /* === Compute set differences ====================================== */ alpar@9: alpar@9: DEBUG3 (("** Computing set differences phase. **\n")) ; alpar@9: alpar@9: /* pivot row is currently dead - it will be revived later. */ alpar@9: alpar@9: DEBUG3 (("Pivot row: ")) ; alpar@9: /* for each column in pivot row */ alpar@9: rp = &A [pivot_row_start] ; alpar@9: rp_end = rp + pivot_row_length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: col = *rp++ ; alpar@9: ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; alpar@9: DEBUG3 (("Col: %d\n", col)) ; alpar@9: alpar@9: /* clear tags used to construct pivot row pattern */ alpar@9: col_thickness = -Col [col].shared1.thickness ; alpar@9: ASSERT (col_thickness > 0) ; alpar@9: Col [col].shared1.thickness = col_thickness ; alpar@9: alpar@9: /* === Remove column from degree list =========================== */ alpar@9: alpar@9: cur_score = Col [col].shared2.score ; alpar@9: prev_col = Col [col].shared3.prev ; alpar@9: next_col = Col [col].shared4.degree_next ; alpar@9: ASSERT (cur_score >= 0) ; alpar@9: ASSERT (cur_score <= n_col) ; alpar@9: ASSERT (cur_score >= EMPTY) ; alpar@9: if (prev_col == EMPTY) alpar@9: { alpar@9: head [cur_score] = next_col ; alpar@9: } alpar@9: else alpar@9: { alpar@9: Col [prev_col].shared4.degree_next = next_col ; alpar@9: } alpar@9: if (next_col != EMPTY) alpar@9: { alpar@9: Col [next_col].shared3.prev = prev_col ; alpar@9: } alpar@9: alpar@9: /* === Scan the column ========================================== */ alpar@9: alpar@9: cp = &A [Col [col].start] ; alpar@9: cp_end = cp + Col [col].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: /* get a row */ alpar@9: row = *cp++ ; alpar@9: row_mark = Row [row].shared2.mark ; alpar@9: /* skip if dead */ alpar@9: if (ROW_IS_MARKED_DEAD (row_mark)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: ASSERT (row != pivot_row) ; alpar@9: set_difference = row_mark - tag_mark ; alpar@9: /* check if the row has been seen yet */ alpar@9: if (set_difference < 0) alpar@9: { alpar@9: ASSERT (Row [row].shared1.degree <= max_deg) ; alpar@9: set_difference = Row [row].shared1.degree ; alpar@9: } alpar@9: /* subtract column thickness from this row's set difference */ alpar@9: set_difference -= col_thickness ; alpar@9: ASSERT (set_difference >= 0) ; alpar@9: /* absorb this row if the set difference becomes zero */ alpar@9: if (set_difference == 0 && aggressive) alpar@9: { alpar@9: DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; alpar@9: KILL_ROW (row) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* save the new mark */ alpar@9: Row [row].shared2.mark = set_difference + tag_mark ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_deg_lists (n_row, n_col, Row, Col, head, alpar@9: min_score, n_col2-k-pivot_row_degree, max_deg) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Add up set differences for each column ======================= */ alpar@9: alpar@9: DEBUG3 (("** Adding set differences phase. **\n")) ; alpar@9: alpar@9: /* for each column in pivot row */ alpar@9: rp = &A [pivot_row_start] ; alpar@9: rp_end = rp + pivot_row_length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: /* get a column */ alpar@9: col = *rp++ ; alpar@9: ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; alpar@9: hash = 0 ; alpar@9: cur_score = 0 ; alpar@9: cp = &A [Col [col].start] ; alpar@9: /* compact the column */ alpar@9: new_cp = cp ; alpar@9: cp_end = cp + Col [col].length ; alpar@9: alpar@9: DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; alpar@9: alpar@9: while (cp < cp_end) alpar@9: { alpar@9: /* get a row */ alpar@9: row = *cp++ ; alpar@9: ASSERT(row >= 0 && row < n_row) ; alpar@9: row_mark = Row [row].shared2.mark ; alpar@9: /* skip if dead */ alpar@9: if (ROW_IS_MARKED_DEAD (row_mark)) alpar@9: { alpar@9: DEBUG4 ((" Row %d, dead\n", row)) ; alpar@9: continue ; alpar@9: } alpar@9: DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark)); alpar@9: ASSERT (row_mark >= tag_mark) ; alpar@9: /* compact the column */ alpar@9: *new_cp++ = row ; alpar@9: /* compute hash function */ alpar@9: hash += row ; alpar@9: /* add set difference */ alpar@9: cur_score += row_mark - tag_mark ; alpar@9: /* integer overflow... */ alpar@9: cur_score = MIN (cur_score, n_col) ; alpar@9: } alpar@9: alpar@9: /* recompute the column's length */ alpar@9: Col [col].length = (Int) (new_cp - &A [Col [col].start]) ; alpar@9: alpar@9: /* === Further mass elimination ================================= */ alpar@9: alpar@9: if (Col [col].length == 0) alpar@9: { alpar@9: DEBUG4 (("further mass elimination. Col: %d\n", col)) ; alpar@9: /* nothing left but the pivot row in this column */ alpar@9: KILL_PRINCIPAL_COL (col) ; alpar@9: pivot_row_degree -= Col [col].shared1.thickness ; alpar@9: ASSERT (pivot_row_degree >= 0) ; alpar@9: /* order it */ alpar@9: Col [col].shared2.order = k ; alpar@9: /* increment order count by column thickness */ alpar@9: k += Col [col].shared1.thickness ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* === Prepare for supercolumn detection ==================== */ alpar@9: alpar@9: DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; alpar@9: alpar@9: /* save score so far */ alpar@9: Col [col].shared2.score = cur_score ; alpar@9: alpar@9: /* add column to hash table, for supercolumn detection */ alpar@9: hash %= n_col + 1 ; alpar@9: alpar@9: DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; alpar@9: ASSERT (((Int) hash) <= n_col) ; alpar@9: alpar@9: head_column = head [hash] ; alpar@9: if (head_column > EMPTY) alpar@9: { alpar@9: /* degree list "hash" is non-empty, use prev (shared3) of */ alpar@9: /* first column in degree list as head of hash bucket */ alpar@9: first_col = Col [head_column].shared3.headhash ; alpar@9: Col [head_column].shared3.headhash = col ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* degree list "hash" is empty, use head as hash bucket */ alpar@9: first_col = - (head_column + 2) ; alpar@9: head [hash] = - (col + 2) ; alpar@9: } alpar@9: Col [col].shared4.hash_next = first_col ; alpar@9: alpar@9: /* save hash function in Col [col].shared3.hash */ alpar@9: Col [col].shared3.hash = (Int) hash ; alpar@9: ASSERT (COL_IS_ALIVE (col)) ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* The approximate external column degree is now computed. */ alpar@9: alpar@9: /* === Supercolumn detection ======================================== */ alpar@9: alpar@9: DEBUG3 (("** Supercolumn detection phase. **\n")) ; alpar@9: alpar@9: detect_super_cols ( alpar@9: alpar@9: #ifndef NDEBUG alpar@9: n_col, Row, alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: Col, A, head, pivot_row_start, pivot_row_length) ; alpar@9: alpar@9: /* === Kill the pivotal column ====================================== */ alpar@9: alpar@9: KILL_PRINCIPAL_COL (pivot_col) ; alpar@9: alpar@9: /* === Clear mark =================================================== */ alpar@9: alpar@9: tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ; alpar@9: alpar@9: #ifndef NDEBUG alpar@9: DEBUG3 (("check3\n")) ; alpar@9: debug_mark (n_row, Row, tag_mark, max_mark) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Finalize the new pivot row, and column scores ================ */ alpar@9: alpar@9: DEBUG3 (("** Finalize scores phase. **\n")) ; alpar@9: alpar@9: /* for each column in pivot row */ alpar@9: rp = &A [pivot_row_start] ; alpar@9: /* compact the pivot row */ alpar@9: new_rp = rp ; alpar@9: rp_end = rp + pivot_row_length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: col = *rp++ ; alpar@9: /* skip dead columns */ alpar@9: if (COL_IS_DEAD (col)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: *new_rp++ = col ; alpar@9: /* add new pivot row to column */ alpar@9: A [Col [col].start + (Col [col].length++)] = pivot_row ; alpar@9: alpar@9: /* retrieve score so far and add on pivot row's degree. */ alpar@9: /* (we wait until here for this in case the pivot */ alpar@9: /* row's degree was reduced due to mass elimination). */ alpar@9: cur_score = Col [col].shared2.score + pivot_row_degree ; alpar@9: alpar@9: /* calculate the max possible score as the number of */ alpar@9: /* external columns minus the 'k' value minus the */ alpar@9: /* columns thickness */ alpar@9: max_score = n_col - k - Col [col].shared1.thickness ; alpar@9: alpar@9: /* make the score the external degree of the union-of-rows */ alpar@9: cur_score -= Col [col].shared1.thickness ; alpar@9: alpar@9: /* make sure score is less or equal than the max score */ alpar@9: cur_score = MIN (cur_score, max_score) ; alpar@9: ASSERT (cur_score >= 0) ; alpar@9: alpar@9: /* store updated score */ alpar@9: Col [col].shared2.score = cur_score ; alpar@9: alpar@9: /* === Place column back in degree list ========================= */ alpar@9: alpar@9: ASSERT (min_score >= 0) ; alpar@9: ASSERT (min_score <= n_col) ; alpar@9: ASSERT (cur_score >= 0) ; alpar@9: ASSERT (cur_score <= n_col) ; alpar@9: ASSERT (head [cur_score] >= EMPTY) ; alpar@9: next_col = head [cur_score] ; alpar@9: Col [col].shared4.degree_next = next_col ; alpar@9: Col [col].shared3.prev = EMPTY ; alpar@9: if (next_col != EMPTY) alpar@9: { alpar@9: Col [next_col].shared3.prev = col ; alpar@9: } alpar@9: head [cur_score] = col ; alpar@9: alpar@9: /* see if this score is less than current min */ alpar@9: min_score = MIN (min_score, cur_score) ; alpar@9: alpar@9: } alpar@9: alpar@9: #ifndef NDEBUG alpar@9: debug_deg_lists (n_row, n_col, Row, Col, head, alpar@9: min_score, n_col2-k, max_deg) ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Resurrect the new pivot row ================================== */ alpar@9: alpar@9: if (pivot_row_degree > 0) alpar@9: { alpar@9: /* update pivot row length to reflect any cols that were killed */ alpar@9: /* during super-col detection and mass elimination */ alpar@9: Row [pivot_row].start = pivot_row_start ; alpar@9: Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ; alpar@9: ASSERT (Row [pivot_row].length > 0) ; alpar@9: Row [pivot_row].shared1.degree = pivot_row_degree ; alpar@9: Row [pivot_row].shared2.mark = 0 ; alpar@9: /* pivot row is no longer dead */ alpar@9: alpar@9: DEBUG1 (("Resurrect Pivot_row %d deg: %d\n", alpar@9: pivot_row, pivot_row_degree)) ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* === All principal columns have now been ordered ====================== */ alpar@9: alpar@9: return (ngarbage) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === order_children ======================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: The find_ordering routine has ordered all of the principal columns (the alpar@9: representatives of the supercolumns). The non-principal columns have not alpar@9: yet been ordered. This routine orders those columns by walking up the alpar@9: parent tree (a column is a child of the column which absorbed it). The alpar@9: final permutation vector is then placed in p [0 ... n_col-1], with p [0] alpar@9: being the first column, and p [n_col-1] being the last. It doesn't look alpar@9: like it at first glance, but be assured that this routine takes time linear alpar@9: in the number of columns. Although not immediately obvious, the time alpar@9: taken by this routine is O (n_col), that is, linear in the number of alpar@9: columns. Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE void order_children alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_col, /* number of columns of A */ alpar@9: Colamd_Col Col [], /* of size n_col+1 */ alpar@9: Int p [] /* p [0 ... n_col-1] is the column permutation*/ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int i ; /* loop counter for all columns */ alpar@9: Int c ; /* column index */ alpar@9: Int parent ; /* index of column's parent */ alpar@9: Int order ; /* column's order */ alpar@9: alpar@9: /* === Order each non-principal column ================================== */ alpar@9: alpar@9: for (i = 0 ; i < n_col ; i++) alpar@9: { alpar@9: /* find an un-ordered non-principal column */ alpar@9: ASSERT (COL_IS_DEAD (i)) ; alpar@9: if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY) alpar@9: { alpar@9: parent = i ; alpar@9: /* once found, find its principal parent */ alpar@9: do alpar@9: { alpar@9: parent = Col [parent].shared1.parent ; alpar@9: } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; alpar@9: alpar@9: /* now, order all un-ordered non-principal columns along path */ alpar@9: /* to this parent. collapse tree at the same time */ alpar@9: c = i ; alpar@9: /* get order of parent */ alpar@9: order = Col [parent].shared2.order ; alpar@9: alpar@9: do alpar@9: { alpar@9: ASSERT (Col [c].shared2.order == EMPTY) ; alpar@9: alpar@9: /* order this column */ alpar@9: Col [c].shared2.order = order++ ; alpar@9: /* collaps tree */ alpar@9: Col [c].shared1.parent = parent ; alpar@9: alpar@9: /* get immediate parent of this column */ alpar@9: c = Col [c].shared1.parent ; alpar@9: alpar@9: /* continue until we hit an ordered column. There are */ alpar@9: /* guarranteed not to be anymore unordered columns */ alpar@9: /* above an ordered column */ alpar@9: } while (Col [c].shared2.order == EMPTY) ; alpar@9: alpar@9: /* re-order the super_col parent to largest order for this group */ alpar@9: Col [parent].shared2.order = order ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Generate the permutation ========================================= */ alpar@9: alpar@9: for (c = 0 ; c < n_col ; c++) alpar@9: { alpar@9: p [Col [c].shared2.order] = c ; alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === detect_super_cols ==================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Detects supercolumns by finding matches between columns in the hash buckets. alpar@9: Check amongst columns in the set A [row_start ... row_start + row_length-1]. alpar@9: The columns under consideration are currently *not* in the degree lists, alpar@9: and have already been placed in the hash buckets. alpar@9: alpar@9: The hash bucket for columns whose hash function is equal to h is stored alpar@9: as follows: alpar@9: alpar@9: if head [h] is >= 0, then head [h] contains a degree list, so: alpar@9: alpar@9: head [h] is the first column in degree bucket h. alpar@9: Col [head [h]].headhash gives the first column in hash bucket h. alpar@9: alpar@9: otherwise, the degree list is empty, and: alpar@9: alpar@9: -(head [h] + 2) is the first column in hash bucket h. alpar@9: alpar@9: For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous alpar@9: column" pointer. Col [c].shared3.hash is used instead as the hash number alpar@9: for that column. The value of Col [c].shared4.hash_next is the next column alpar@9: in the same hash bucket. alpar@9: alpar@9: Assuming no, or "few" hash collisions, the time taken by this routine is alpar@9: linear in the sum of the sizes (lengths) of each column whose score has alpar@9: just been computed in the approximate degree computation. alpar@9: Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE void detect_super_cols alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: /* these two parameters are only needed when debugging is enabled: */ alpar@9: Int n_col, /* number of columns of A */ alpar@9: Colamd_Row Row [], /* of size n_row+1 */ alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: Colamd_Col Col [], /* of size n_col+1 */ alpar@9: Int A [], /* row indices of A */ alpar@9: Int head [], /* head of degree lists and hash buckets */ alpar@9: Int row_start, /* pointer to set of columns to check */ alpar@9: Int row_length /* number of columns to check */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int hash ; /* hash value for a column */ alpar@9: Int *rp ; /* pointer to a row */ alpar@9: Int c ; /* a column index */ alpar@9: Int super_c ; /* column index of the column to absorb into */ alpar@9: Int *cp1 ; /* column pointer for column super_c */ alpar@9: Int *cp2 ; /* column pointer for column c */ alpar@9: Int length ; /* length of column super_c */ alpar@9: Int prev_c ; /* column preceding c in hash bucket */ alpar@9: Int i ; /* loop counter */ alpar@9: Int *rp_end ; /* pointer to the end of the row */ alpar@9: Int col ; /* a column index in the row to check */ alpar@9: Int head_column ; /* first column in hash bucket or degree list */ alpar@9: Int first_col ; /* first column in hash bucket */ alpar@9: alpar@9: /* === Consider each column in the row ================================== */ alpar@9: alpar@9: rp = &A [row_start] ; alpar@9: rp_end = rp + row_length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: col = *rp++ ; alpar@9: if (COL_IS_DEAD (col)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: alpar@9: /* get hash number for this column */ alpar@9: hash = Col [col].shared3.hash ; alpar@9: ASSERT (hash <= n_col) ; alpar@9: alpar@9: /* === Get the first column in this hash bucket ===================== */ alpar@9: alpar@9: head_column = head [hash] ; alpar@9: if (head_column > EMPTY) alpar@9: { alpar@9: first_col = Col [head_column].shared3.headhash ; alpar@9: } alpar@9: else alpar@9: { alpar@9: first_col = - (head_column + 2) ; alpar@9: } alpar@9: alpar@9: /* === Consider each column in the hash bucket ====================== */ alpar@9: alpar@9: for (super_c = first_col ; super_c != EMPTY ; alpar@9: super_c = Col [super_c].shared4.hash_next) alpar@9: { alpar@9: ASSERT (COL_IS_ALIVE (super_c)) ; alpar@9: ASSERT (Col [super_c].shared3.hash == hash) ; alpar@9: length = Col [super_c].length ; alpar@9: alpar@9: /* prev_c is the column preceding column c in the hash bucket */ alpar@9: prev_c = super_c ; alpar@9: alpar@9: /* === Compare super_c with all columns after it ================ */ alpar@9: alpar@9: for (c = Col [super_c].shared4.hash_next ; alpar@9: c != EMPTY ; c = Col [c].shared4.hash_next) alpar@9: { alpar@9: ASSERT (c != super_c) ; alpar@9: ASSERT (COL_IS_ALIVE (c)) ; alpar@9: ASSERT (Col [c].shared3.hash == hash) ; alpar@9: alpar@9: /* not identical if lengths or scores are different */ alpar@9: if (Col [c].length != length || alpar@9: Col [c].shared2.score != Col [super_c].shared2.score) alpar@9: { alpar@9: prev_c = c ; alpar@9: continue ; alpar@9: } alpar@9: alpar@9: /* compare the two columns */ alpar@9: cp1 = &A [Col [super_c].start] ; alpar@9: cp2 = &A [Col [c].start] ; alpar@9: alpar@9: for (i = 0 ; i < length ; i++) alpar@9: { alpar@9: /* the columns are "clean" (no dead rows) */ alpar@9: ASSERT (ROW_IS_ALIVE (*cp1)) ; alpar@9: ASSERT (ROW_IS_ALIVE (*cp2)) ; alpar@9: /* row indices will same order for both supercols, */ alpar@9: /* no gather scatter nessasary */ alpar@9: if (*cp1++ != *cp2++) alpar@9: { alpar@9: break ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* the two columns are different if the for-loop "broke" */ alpar@9: if (i != length) alpar@9: { alpar@9: prev_c = c ; alpar@9: continue ; alpar@9: } alpar@9: alpar@9: /* === Got it! two columns are identical =================== */ alpar@9: alpar@9: ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; alpar@9: alpar@9: Col [super_c].shared1.thickness += Col [c].shared1.thickness ; alpar@9: Col [c].shared1.parent = super_c ; alpar@9: KILL_NON_PRINCIPAL_COL (c) ; alpar@9: /* order c later, in order_children() */ alpar@9: Col [c].shared2.order = EMPTY ; alpar@9: /* remove c from hash bucket */ alpar@9: Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Empty this hash bucket ======================================= */ alpar@9: alpar@9: if (head_column > EMPTY) alpar@9: { alpar@9: /* corresponding degree list "hash" is not empty */ alpar@9: Col [head_column].shared3.headhash = EMPTY ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* corresponding degree list "hash" is empty */ alpar@9: head [hash] = EMPTY ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === garbage_collection =================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Defragments and compacts columns and rows in the workspace A. Used when alpar@9: all avaliable memory has been used while performing row merging. Returns alpar@9: the index of the first free position in A, after garbage collection. The alpar@9: time taken by this routine is linear is the size of the array A, which is alpar@9: itself linear in the number of nonzeros in the input matrix. alpar@9: Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE Int garbage_collection /* returns the new value of pfree */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, /* number of rows */ alpar@9: Int n_col, /* number of columns */ alpar@9: Colamd_Row Row [], /* row info */ alpar@9: Colamd_Col Col [], /* column info */ alpar@9: Int A [], /* A [0 ... Alen-1] holds the matrix */ alpar@9: Int *pfree /* &A [0] ... pfree is in use */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int *psrc ; /* source pointer */ alpar@9: Int *pdest ; /* destination pointer */ alpar@9: Int j ; /* counter */ alpar@9: Int r ; /* a row index */ alpar@9: Int c ; /* a column index */ alpar@9: Int length ; /* length of a row or column */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: Int debug_rows ; alpar@9: DEBUG2 (("Defrag..\n")) ; alpar@9: for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ; alpar@9: debug_rows = 0 ; alpar@9: #endif /* NDEBUG */ alpar@9: alpar@9: /* === Defragment the columns =========================================== */ alpar@9: alpar@9: pdest = &A[0] ; alpar@9: for (c = 0 ; c < n_col ; c++) alpar@9: { alpar@9: if (COL_IS_ALIVE (c)) alpar@9: { alpar@9: psrc = &A [Col [c].start] ; alpar@9: alpar@9: /* move and compact the column */ alpar@9: ASSERT (pdest <= psrc) ; alpar@9: Col [c].start = (Int) (pdest - &A [0]) ; alpar@9: length = Col [c].length ; alpar@9: for (j = 0 ; j < length ; j++) alpar@9: { alpar@9: r = *psrc++ ; alpar@9: if (ROW_IS_ALIVE (r)) alpar@9: { alpar@9: *pdest++ = r ; alpar@9: } alpar@9: } alpar@9: Col [c].length = (Int) (pdest - &A [Col [c].start]) ; alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Prepare to defragment the rows =================================== */ alpar@9: alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: if (ROW_IS_DEAD (r) || (Row [r].length == 0)) alpar@9: { alpar@9: /* This row is already dead, or is of zero length. Cannot compact alpar@9: * a row of zero length, so kill it. NOTE: in the current version, alpar@9: * there are no zero-length live rows. Kill the row (for the first alpar@9: * time, or again) just to be safe. */ alpar@9: KILL_ROW (r) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: /* save first column index in Row [r].shared2.first_column */ alpar@9: psrc = &A [Row [r].start] ; alpar@9: Row [r].shared2.first_column = *psrc ; alpar@9: ASSERT (ROW_IS_ALIVE (r)) ; alpar@9: /* flag the start of the row with the one's complement of row */ alpar@9: *psrc = ONES_COMPLEMENT (r) ; alpar@9: #ifndef NDEBUG alpar@9: debug_rows++ ; alpar@9: #endif /* NDEBUG */ alpar@9: } alpar@9: } alpar@9: alpar@9: /* === Defragment the rows ============================================== */ alpar@9: alpar@9: psrc = pdest ; alpar@9: while (psrc < pfree) alpar@9: { alpar@9: /* find a negative number ... the start of a row */ alpar@9: if (*psrc++ < 0) alpar@9: { alpar@9: psrc-- ; alpar@9: /* get the row index */ alpar@9: r = ONES_COMPLEMENT (*psrc) ; alpar@9: ASSERT (r >= 0 && r < n_row) ; alpar@9: /* restore first column index */ alpar@9: *psrc = Row [r].shared2.first_column ; alpar@9: ASSERT (ROW_IS_ALIVE (r)) ; alpar@9: ASSERT (Row [r].length > 0) ; alpar@9: /* move and compact the row */ alpar@9: ASSERT (pdest <= psrc) ; alpar@9: Row [r].start = (Int) (pdest - &A [0]) ; alpar@9: length = Row [r].length ; alpar@9: for (j = 0 ; j < length ; j++) alpar@9: { alpar@9: c = *psrc++ ; alpar@9: if (COL_IS_ALIVE (c)) alpar@9: { alpar@9: *pdest++ = c ; alpar@9: } alpar@9: } alpar@9: Row [r].length = (Int) (pdest - &A [Row [r].start]) ; alpar@9: ASSERT (Row [r].length > 0) ; alpar@9: #ifndef NDEBUG alpar@9: debug_rows-- ; alpar@9: #endif /* NDEBUG */ alpar@9: } alpar@9: } alpar@9: /* ensure we found all the rows */ alpar@9: ASSERT (debug_rows == 0) ; alpar@9: alpar@9: /* === Return the new value of pfree ==================================== */ alpar@9: alpar@9: return ((Int) (pdest - &A [0])) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === clear_mark =========================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Clears the Row [].shared2.mark array, and returns the new tag_mark. alpar@9: Return value is the new tag_mark. Not user-callable. alpar@9: */ alpar@9: alpar@9: PRIVATE Int clear_mark /* return the new value for tag_mark */ alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int tag_mark, /* new value of tag_mark */ alpar@9: Int max_mark, /* max allowed value of tag_mark */ alpar@9: alpar@9: Int n_row, /* number of rows in A */ alpar@9: Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int r ; alpar@9: alpar@9: if (tag_mark <= 0 || tag_mark >= max_mark) alpar@9: { alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: if (ROW_IS_ALIVE (r)) alpar@9: { alpar@9: Row [r].shared2.mark = 0 ; alpar@9: } alpar@9: } alpar@9: tag_mark = 1 ; alpar@9: } alpar@9: alpar@9: return (tag_mark) ; alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === print_report ========================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: PRIVATE void print_report alpar@9: ( alpar@9: char *method, alpar@9: Int stats [COLAMD_STATS] alpar@9: ) alpar@9: { alpar@9: alpar@9: Int i1, i2, i3 ; alpar@9: alpar@9: PRINTF (("\n%s version %d.%d, %s: ", method, alpar@9: COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ; alpar@9: alpar@9: if (!stats) alpar@9: { alpar@9: PRINTF (("No statistics available.\n")) ; alpar@9: return ; alpar@9: } alpar@9: alpar@9: i1 = stats [COLAMD_INFO1] ; alpar@9: i2 = stats [COLAMD_INFO2] ; alpar@9: i3 = stats [COLAMD_INFO3] ; alpar@9: alpar@9: if (stats [COLAMD_STATUS] >= 0) alpar@9: { alpar@9: PRINTF (("OK. ")) ; alpar@9: } alpar@9: else alpar@9: { alpar@9: PRINTF (("ERROR. ")) ; alpar@9: } alpar@9: alpar@9: switch (stats [COLAMD_STATUS]) alpar@9: { alpar@9: alpar@9: case COLAMD_OK_BUT_JUMBLED: alpar@9: alpar@9: PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ; alpar@9: alpar@9: PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n", alpar@9: method, i3)) ; alpar@9: alpar@9: PRINTF(("%s: last seen duplicate or out-of-order row index: %d\n", alpar@9: method, INDEX (i2))) ; alpar@9: alpar@9: PRINTF(("%s: last seen in column: %d", alpar@9: method, INDEX (i1))) ; alpar@9: alpar@9: /* no break - fall through to next case instead */ alpar@9: alpar@9: case COLAMD_OK: alpar@9: alpar@9: PRINTF(("\n")) ; alpar@9: alpar@9: PRINTF(("%s: number of dense or empty rows ignored: %d\n", alpar@9: method, stats [COLAMD_DENSE_ROW])) ; alpar@9: alpar@9: PRINTF(("%s: number of dense or empty columns ignored: %d\n", alpar@9: method, stats [COLAMD_DENSE_COL])) ; alpar@9: alpar@9: PRINTF(("%s: number of garbage collections performed: %d\n", alpar@9: method, stats [COLAMD_DEFRAG_COUNT])) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_A_not_present: alpar@9: alpar@9: PRINTF(("Array A (row indices of matrix) not present.\n")) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_p_not_present: alpar@9: alpar@9: PRINTF(("Array p (column pointers for matrix) not present.\n")) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_nrow_negative: alpar@9: alpar@9: PRINTF(("Invalid number of rows (%d).\n", i1)) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_ncol_negative: alpar@9: alpar@9: PRINTF(("Invalid number of columns (%d).\n", i1)) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_nnz_negative: alpar@9: alpar@9: PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_p0_nonzero: alpar@9: alpar@9: PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1)); alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_A_too_small: alpar@9: alpar@9: PRINTF(("Array A too small.\n")) ; alpar@9: PRINTF((" Need Alen >= %d, but given only Alen = %d.\n", alpar@9: i1, i2)) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_col_length_negative: alpar@9: alpar@9: PRINTF alpar@9: (("Column %d has a negative number of nonzero entries (%d).\n", alpar@9: INDEX (i1), i2)) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_row_index_out_of_bounds: alpar@9: alpar@9: PRINTF alpar@9: (("Row index (row %d) out of bounds (%d to %d) in column %d.\n", alpar@9: INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ; alpar@9: break ; alpar@9: alpar@9: case COLAMD_ERROR_out_of_memory: alpar@9: alpar@9: PRINTF(("Out of memory.\n")) ; alpar@9: break ; alpar@9: alpar@9: /* v2.4: internal-error case deleted */ alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === colamd debugging routines ============================================ */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* When debugging is disabled, the remainder of this file is ignored. */ alpar@9: alpar@9: #ifndef NDEBUG alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === debug_structures ===================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: At this point, all empty rows and columns are dead. All live columns alpar@9: are "clean" (containing no dead rows) and simplicial (no supercolumns alpar@9: yet). Rows may contain dead columns, but all live rows contain at alpar@9: least one live column. alpar@9: */ alpar@9: alpar@9: PRIVATE void debug_structures alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [], alpar@9: Int n_col2 alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int i ; alpar@9: Int c ; alpar@9: Int *cp ; alpar@9: Int *cp_end ; alpar@9: Int len ; alpar@9: Int score ; alpar@9: Int r ; alpar@9: Int *rp ; alpar@9: Int *rp_end ; alpar@9: Int deg ; alpar@9: alpar@9: /* === Check A, Row, and Col ============================================ */ alpar@9: alpar@9: for (c = 0 ; c < n_col ; c++) alpar@9: { alpar@9: if (COL_IS_ALIVE (c)) alpar@9: { alpar@9: len = Col [c].length ; alpar@9: score = Col [c].shared2.score ; alpar@9: DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; alpar@9: ASSERT (len > 0) ; alpar@9: ASSERT (score >= 0) ; alpar@9: ASSERT (Col [c].shared1.thickness == 1) ; alpar@9: cp = &A [Col [c].start] ; alpar@9: cp_end = cp + len ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: r = *cp++ ; alpar@9: ASSERT (ROW_IS_ALIVE (r)) ; alpar@9: } alpar@9: } alpar@9: else alpar@9: { alpar@9: i = Col [c].shared2.order ; alpar@9: ASSERT (i >= n_col2 && i < n_col) ; alpar@9: } alpar@9: } alpar@9: alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: if (ROW_IS_ALIVE (r)) alpar@9: { alpar@9: i = 0 ; alpar@9: len = Row [r].length ; alpar@9: deg = Row [r].shared1.degree ; alpar@9: ASSERT (len > 0) ; alpar@9: ASSERT (deg > 0) ; alpar@9: rp = &A [Row [r].start] ; alpar@9: rp_end = rp + len ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: c = *rp++ ; alpar@9: if (COL_IS_ALIVE (c)) alpar@9: { alpar@9: i++ ; alpar@9: } alpar@9: } alpar@9: ASSERT (i > 0) ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === debug_deg_lists ====================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Prints the contents of the degree lists. Counts the number of columns alpar@9: in the degree list and compares it to the total it should have. Also alpar@9: checks the row degrees. alpar@9: */ alpar@9: alpar@9: PRIVATE void debug_deg_lists alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int head [], alpar@9: Int min_score, alpar@9: Int should, alpar@9: Int max_deg alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int deg ; alpar@9: Int col ; alpar@9: Int have ; alpar@9: Int row ; alpar@9: alpar@9: /* === Check the degree lists =========================================== */ alpar@9: alpar@9: if (n_col > 10000 && colamd_debug <= 0) alpar@9: { alpar@9: return ; alpar@9: } alpar@9: have = 0 ; alpar@9: DEBUG4 (("Degree lists: %d\n", min_score)) ; alpar@9: for (deg = 0 ; deg <= n_col ; deg++) alpar@9: { alpar@9: col = head [deg] ; alpar@9: if (col == EMPTY) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: DEBUG4 (("%d:", deg)) ; alpar@9: while (col != EMPTY) alpar@9: { alpar@9: DEBUG4 ((" %d", col)) ; alpar@9: have += Col [col].shared1.thickness ; alpar@9: ASSERT (COL_IS_ALIVE (col)) ; alpar@9: col = Col [col].shared4.degree_next ; alpar@9: } alpar@9: DEBUG4 (("\n")) ; alpar@9: } alpar@9: DEBUG4 (("should %d have %d\n", should, have)) ; alpar@9: ASSERT (should == have) ; alpar@9: alpar@9: /* === Check the row degrees ============================================ */ alpar@9: alpar@9: if (n_row > 10000 && colamd_debug <= 0) alpar@9: { alpar@9: return ; alpar@9: } alpar@9: for (row = 0 ; row < n_row ; row++) alpar@9: { alpar@9: if (ROW_IS_ALIVE (row)) alpar@9: { alpar@9: ASSERT (Row [row].shared1.degree <= max_deg) ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === debug_mark =========================================================== */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Ensures that the tag_mark is less that the maximum and also ensures that alpar@9: each entry in the mark array is less than the tag mark. alpar@9: */ alpar@9: alpar@9: PRIVATE void debug_mark alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, alpar@9: Colamd_Row Row [], alpar@9: Int tag_mark, alpar@9: Int max_mark alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int r ; alpar@9: alpar@9: /* === Check the Row marks ============================================== */ alpar@9: alpar@9: ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; alpar@9: if (n_row > 10000 && colamd_debug <= 0) alpar@9: { alpar@9: return ; alpar@9: } alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: ASSERT (Row [r].shared2.mark < tag_mark) ; alpar@9: } alpar@9: } alpar@9: alpar@9: alpar@9: /* ========================================================================== */ alpar@9: /* === debug_matrix ========================================================= */ alpar@9: /* ========================================================================== */ alpar@9: alpar@9: /* alpar@9: Prints out the contents of the columns and the rows. alpar@9: */ alpar@9: alpar@9: PRIVATE void debug_matrix alpar@9: ( alpar@9: /* === Parameters ======================================================= */ alpar@9: alpar@9: Int n_row, alpar@9: Int n_col, alpar@9: Colamd_Row Row [], alpar@9: Colamd_Col Col [], alpar@9: Int A [] alpar@9: ) alpar@9: { alpar@9: /* === Local variables ================================================== */ alpar@9: alpar@9: Int r ; alpar@9: Int c ; alpar@9: Int *rp ; alpar@9: Int *rp_end ; alpar@9: Int *cp ; alpar@9: Int *cp_end ; alpar@9: alpar@9: /* === Dump the rows and columns of the matrix ========================== */ alpar@9: alpar@9: if (colamd_debug < 3) alpar@9: { alpar@9: return ; alpar@9: } alpar@9: DEBUG3 (("DUMP MATRIX:\n")) ; alpar@9: for (r = 0 ; r < n_row ; r++) alpar@9: { alpar@9: DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ; alpar@9: if (ROW_IS_DEAD (r)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: DEBUG3 (("start %d length %d degree %d\n", alpar@9: Row [r].start, Row [r].length, Row [r].shared1.degree)) ; alpar@9: rp = &A [Row [r].start] ; alpar@9: rp_end = rp + Row [r].length ; alpar@9: while (rp < rp_end) alpar@9: { alpar@9: c = *rp++ ; alpar@9: DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ; alpar@9: } alpar@9: } alpar@9: alpar@9: for (c = 0 ; c < n_col ; c++) alpar@9: { alpar@9: DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ; alpar@9: if (COL_IS_DEAD (c)) alpar@9: { alpar@9: continue ; alpar@9: } alpar@9: DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", alpar@9: Col [c].start, Col [c].length, alpar@9: Col [c].shared1.thickness, Col [c].shared2.score)) ; alpar@9: cp = &A [Col [c].start] ; alpar@9: cp_end = cp + Col [c].length ; alpar@9: while (cp < cp_end) alpar@9: { alpar@9: r = *cp++ ; alpar@9: DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ; alpar@9: } alpar@9: } alpar@9: } alpar@9: alpar@9: PRIVATE void colamd_get_debug alpar@9: ( alpar@9: char *method alpar@9: ) alpar@9: { alpar@9: FILE *f ; alpar@9: colamd_debug = 0 ; /* no debug printing */ alpar@9: f = fopen ("debug", "r") ; alpar@9: if (f == (FILE *) NULL) alpar@9: { alpar@9: colamd_debug = 0 ; alpar@9: } alpar@9: else alpar@9: { alpar@9: fscanf (f, "%d", &colamd_debug) ; alpar@9: fclose (f) ; alpar@9: } alpar@9: DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", alpar@9: method, colamd_debug)) ; alpar@9: } alpar@9: alpar@9: #endif /* NDEBUG */