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