1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/colamd/colamd.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,3622 @@
1.4 +/* ========================================================================== */
1.5 +/* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
1.6 +/* ========================================================================== */
1.7 +
1.8 +/* COLAMD / SYMAMD
1.9 +
1.10 + colamd: an approximate minimum degree column ordering algorithm,
1.11 + for LU factorization of symmetric or unsymmetric matrices,
1.12 + QR factorization, least squares, interior point methods for
1.13 + linear programming problems, and other related problems.
1.14 +
1.15 + symamd: an approximate minimum degree ordering algorithm for Cholesky
1.16 + factorization of symmetric matrices.
1.17 +
1.18 + Purpose:
1.19 +
1.20 + Colamd computes a permutation Q such that the Cholesky factorization of
1.21 + (AQ)'(AQ) has less fill-in and requires fewer floating point operations
1.22 + than A'A. This also provides a good ordering for sparse partial
1.23 + pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
1.24 + factorization, and P is computed during numerical factorization via
1.25 + conventional partial pivoting with row interchanges. Colamd is the
1.26 + column ordering method used in SuperLU, part of the ScaLAPACK library.
1.27 + It is also available as built-in function in MATLAB Version 6,
1.28 + available from MathWorks, Inc. (http://www.mathworks.com). This
1.29 + routine can be used in place of colmmd in MATLAB.
1.30 +
1.31 + Symamd computes a permutation P of a symmetric matrix A such that the
1.32 + Cholesky factorization of PAP' has less fill-in and requires fewer
1.33 + floating point operations than A. Symamd constructs a matrix M such
1.34 + that M'M has the same nonzero pattern of A, and then orders the columns
1.35 + of M using colmmd. The column ordering of M is then returned as the
1.36 + row and column ordering P of A.
1.37 +
1.38 + Authors:
1.39 +
1.40 + The authors of the code itself are Stefan I. Larimore and Timothy A.
1.41 + Davis (davis at cise.ufl.edu), University of Florida. The algorithm was
1.42 + developed in collaboration with John Gilbert, Xerox PARC, and Esmond
1.43 + Ng, Oak Ridge National Laboratory.
1.44 +
1.45 + Acknowledgements:
1.46 +
1.47 + This work was supported by the National Science Foundation, under
1.48 + grants DMS-9504974 and DMS-9803599.
1.49 +
1.50 + Copyright and License:
1.51 +
1.52 + Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
1.53 + COLAMD is also available under alternate licenses, contact T. Davis
1.54 + for details.
1.55 +
1.56 + This library is free software; you can redistribute it and/or
1.57 + modify it under the terms of the GNU Lesser General Public
1.58 + License as published by the Free Software Foundation; either
1.59 + version 2.1 of the License, or (at your option) any later version.
1.60 +
1.61 + This library is distributed in the hope that it will be useful,
1.62 + but WITHOUT ANY WARRANTY; without even the implied warranty of
1.63 + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1.64 + Lesser General Public License for more details.
1.65 +
1.66 + You should have received a copy of the GNU Lesser General Public
1.67 + License along with this library; if not, write to the Free Software
1.68 + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
1.69 + USA
1.70 +
1.71 + Permission is hereby granted to use or copy this program under the
1.72 + terms of the GNU LGPL, provided that the Copyright, this License,
1.73 + and the Availability of the original version is retained on all copies.
1.74 + User documentation of any code that uses this code or any modified
1.75 + version of this code must cite the Copyright, this License, the
1.76 + Availability note, and "Used by permission." Permission to modify
1.77 + the code and to distribute modified code is granted, provided the
1.78 + Copyright, this License, and the Availability note are retained,
1.79 + and a notice that the code was modified is included.
1.80 +
1.81 + Availability:
1.82 +
1.83 + The colamd/symamd library is available at
1.84 +
1.85 + http://www.cise.ufl.edu/research/sparse/colamd/
1.86 +
1.87 + This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
1.88 + file. It requires the colamd.h file. It is required by the colamdmex.c
1.89 + and symamdmex.c files, for the MATLAB interface to colamd and symamd.
1.90 + Appears as ACM Algorithm 836.
1.91 +
1.92 + See the ChangeLog file for changes since Version 1.0.
1.93 +
1.94 + References:
1.95 +
1.96 + T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
1.97 + minimum degree ordering algorithm, ACM Transactions on Mathematical
1.98 + Software, vol. 30, no. 3., pp. 353-376, 2004.
1.99 +
1.100 + T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
1.101 + an approximate column minimum degree ordering algorithm, ACM
1.102 + Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
1.103 + 2004.
1.104 +
1.105 +*/
1.106 +
1.107 +/* ========================================================================== */
1.108 +/* === Description of user-callable routines ================================ */
1.109 +/* ========================================================================== */
1.110 +
1.111 +/* COLAMD includes both int and UF_long versions of all its routines. The
1.112 + * description below is for the int version. For UF_long, all int arguments
1.113 + * become UF_long. UF_long is normally defined as long, except for WIN64.
1.114 +
1.115 + ----------------------------------------------------------------------------
1.116 + colamd_recommended:
1.117 + ----------------------------------------------------------------------------
1.118 +
1.119 + C syntax:
1.120 +
1.121 + #include "colamd.h"
1.122 + size_t colamd_recommended (int nnz, int n_row, int n_col) ;
1.123 + size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
1.124 + UF_long n_col) ;
1.125 +
1.126 + Purpose:
1.127 +
1.128 + Returns recommended value of Alen for use by colamd. Returns 0
1.129 + if any input argument is negative. The use of this routine
1.130 + is optional. Not needed for symamd, which dynamically allocates
1.131 + its own memory.
1.132 +
1.133 + Note that in v2.4 and earlier, these routines returned int or long.
1.134 + They now return a value of type size_t.
1.135 +
1.136 + Arguments (all input arguments):
1.137 +
1.138 + int nnz ; Number of nonzeros in the matrix A. This must
1.139 + be the same value as p [n_col] in the call to
1.140 + colamd - otherwise you will get a wrong value
1.141 + of the recommended memory to use.
1.142 +
1.143 + int n_row ; Number of rows in the matrix A.
1.144 +
1.145 + int n_col ; Number of columns in the matrix A.
1.146 +
1.147 + ----------------------------------------------------------------------------
1.148 + colamd_set_defaults:
1.149 + ----------------------------------------------------------------------------
1.150 +
1.151 + C syntax:
1.152 +
1.153 + #include "colamd.h"
1.154 + colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
1.155 + colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
1.156 +
1.157 + Purpose:
1.158 +
1.159 + Sets the default parameters. The use of this routine is optional.
1.160 +
1.161 + Arguments:
1.162 +
1.163 + double knobs [COLAMD_KNOBS] ; Output only.
1.164 +
1.165 + NOTE: the meaning of the dense row/col knobs has changed in v2.4
1.166 +
1.167 + knobs [0] and knobs [1] control dense row and col detection:
1.168 +
1.169 + Colamd: rows with more than
1.170 + max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
1.171 + entries are removed prior to ordering. Columns with more than
1.172 + max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
1.173 + entries are removed prior to
1.174 + ordering, and placed last in the output column ordering.
1.175 +
1.176 + Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
1.177 + Rows and columns with more than
1.178 + max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
1.179 + entries are removed prior to ordering, and placed last in the
1.180 + output ordering.
1.181 +
1.182 + COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
1.183 + respectively, in colamd.h. Default values of these two knobs
1.184 + are both 10. Currently, only knobs [0] and knobs [1] are
1.185 + used, but future versions may use more knobs. If so, they will
1.186 + be properly set to their defaults by the future version of
1.187 + colamd_set_defaults, so that the code that calls colamd will
1.188 + not need to change, assuming that you either use
1.189 + colamd_set_defaults, or pass a (double *) NULL pointer as the
1.190 + knobs array to colamd or symamd.
1.191 +
1.192 + knobs [2]: aggressive absorption
1.193 +
1.194 + knobs [COLAMD_AGGRESSIVE] controls whether or not to do
1.195 + aggressive absorption during the ordering. Default is TRUE.
1.196 +
1.197 +
1.198 + ----------------------------------------------------------------------------
1.199 + colamd:
1.200 + ----------------------------------------------------------------------------
1.201 +
1.202 + C syntax:
1.203 +
1.204 + #include "colamd.h"
1.205 + int colamd (int n_row, int n_col, int Alen, int *A, int *p,
1.206 + double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
1.207 + UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
1.208 + UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
1.209 + UF_long stats [COLAMD_STATS]) ;
1.210 +
1.211 + Purpose:
1.212 +
1.213 + Computes a column ordering (Q) of A such that P(AQ)=LU or
1.214 + (AQ)'AQ=LL' have less fill-in and require fewer floating point
1.215 + operations than factorizing the unpermuted matrix A or A'A,
1.216 + respectively.
1.217 +
1.218 + Returns:
1.219 +
1.220 + TRUE (1) if successful, FALSE (0) otherwise.
1.221 +
1.222 + Arguments:
1.223 +
1.224 + int n_row ; Input argument.
1.225 +
1.226 + Number of rows in the matrix A.
1.227 + Restriction: n_row >= 0.
1.228 + Colamd returns FALSE if n_row is negative.
1.229 +
1.230 + int n_col ; Input argument.
1.231 +
1.232 + Number of columns in the matrix A.
1.233 + Restriction: n_col >= 0.
1.234 + Colamd returns FALSE if n_col is negative.
1.235 +
1.236 + int Alen ; Input argument.
1.237 +
1.238 + Restriction (see note):
1.239 + Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
1.240 + Colamd returns FALSE if these conditions are not met.
1.241 +
1.242 + Note: this restriction makes an modest assumption regarding
1.243 + the size of the two typedef's structures in colamd.h.
1.244 + We do, however, guarantee that
1.245 +
1.246 + Alen >= colamd_recommended (nnz, n_row, n_col)
1.247 +
1.248 + will be sufficient. Note: the macro version does not check
1.249 + for integer overflow, and thus is not recommended. Use
1.250 + the colamd_recommended routine instead.
1.251 +
1.252 + int A [Alen] ; Input argument, undefined on output.
1.253 +
1.254 + A is an integer array of size Alen. Alen must be at least as
1.255 + large as the bare minimum value given above, but this is very
1.256 + low, and can result in excessive run time. For best
1.257 + performance, we recommend that Alen be greater than or equal to
1.258 + colamd_recommended (nnz, n_row, n_col), which adds
1.259 + nnz/5 to the bare minimum value given above.
1.260 +
1.261 + On input, the row indices of the entries in column c of the
1.262 + matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices
1.263 + in a given column c need not be in ascending order, and
1.264 + duplicate row indices may be be present. However, colamd will
1.265 + work a little faster if both of these conditions are met
1.266 + (Colamd puts the matrix into this format, if it finds that the
1.267 + the conditions are not met).
1.268 +
1.269 + The matrix is 0-based. That is, rows are in the range 0 to
1.270 + n_row-1, and columns are in the range 0 to n_col-1. Colamd
1.271 + returns FALSE if any row index is out of range.
1.272 +
1.273 + The contents of A are modified during ordering, and are
1.274 + undefined on output.
1.275 +
1.276 + int p [n_col+1] ; Both input and output argument.
1.277 +
1.278 + p is an integer array of size n_col+1. On input, it holds the
1.279 + "pointers" for the column form of the matrix A. Column c of
1.280 + the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
1.281 + entry, p [0], must be zero, and p [c] <= p [c+1] must hold
1.282 + for all c in the range 0 to n_col-1. The value p [n_col] is
1.283 + thus the total number of entries in the pattern of the matrix A.
1.284 + Colamd returns FALSE if these conditions are not met.
1.285 +
1.286 + On output, if colamd returns TRUE, the array p holds the column
1.287 + permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
1.288 + the first column index in the new ordering, and p [n_col-1] is
1.289 + the last. That is, p [k] = j means that column j of A is the
1.290 + kth pivot column, in AQ, where k is in the range 0 to n_col-1
1.291 + (p [0] = j means that column j of A is the first column in AQ).
1.292 +
1.293 + If colamd returns FALSE, then no permutation is returned, and
1.294 + p is undefined on output.
1.295 +
1.296 + double knobs [COLAMD_KNOBS] ; Input argument.
1.297 +
1.298 + See colamd_set_defaults for a description.
1.299 +
1.300 + int stats [COLAMD_STATS] ; Output argument.
1.301 +
1.302 + Statistics on the ordering, and error status.
1.303 + See colamd.h for related definitions.
1.304 + Colamd returns FALSE if stats is not present.
1.305 +
1.306 + stats [0]: number of dense or empty rows ignored.
1.307 +
1.308 + stats [1]: number of dense or empty columns ignored (and
1.309 + ordered last in the output permutation p)
1.310 + Note that a row can become "empty" if it
1.311 + contains only "dense" and/or "empty" columns,
1.312 + and similarly a column can become "empty" if it
1.313 + only contains "dense" and/or "empty" rows.
1.314 +
1.315 + stats [2]: number of garbage collections performed.
1.316 + This can be excessively high if Alen is close
1.317 + to the minimum required value.
1.318 +
1.319 + stats [3]: status code. < 0 is an error code.
1.320 + > 1 is a warning or notice.
1.321 +
1.322 + 0 OK. Each column of the input matrix contained
1.323 + row indices in increasing order, with no
1.324 + duplicates.
1.325 +
1.326 + 1 OK, but columns of input matrix were jumbled
1.327 + (unsorted columns or duplicate entries). Colamd
1.328 + had to do some extra work to sort the matrix
1.329 + first and remove duplicate entries, but it
1.330 + still was able to return a valid permutation
1.331 + (return value of colamd was TRUE).
1.332 +
1.333 + stats [4]: highest numbered column that
1.334 + is unsorted or has duplicate
1.335 + entries.
1.336 + stats [5]: last seen duplicate or
1.337 + unsorted row index.
1.338 + stats [6]: number of duplicate or
1.339 + unsorted row indices.
1.340 +
1.341 + -1 A is a null pointer
1.342 +
1.343 + -2 p is a null pointer
1.344 +
1.345 + -3 n_row is negative
1.346 +
1.347 + stats [4]: n_row
1.348 +
1.349 + -4 n_col is negative
1.350 +
1.351 + stats [4]: n_col
1.352 +
1.353 + -5 number of nonzeros in matrix is negative
1.354 +
1.355 + stats [4]: number of nonzeros, p [n_col]
1.356 +
1.357 + -6 p [0] is nonzero
1.358 +
1.359 + stats [4]: p [0]
1.360 +
1.361 + -7 A is too small
1.362 +
1.363 + stats [4]: required size
1.364 + stats [5]: actual size (Alen)
1.365 +
1.366 + -8 a column has a negative number of entries
1.367 +
1.368 + stats [4]: column with < 0 entries
1.369 + stats [5]: number of entries in col
1.370 +
1.371 + -9 a row index is out of bounds
1.372 +
1.373 + stats [4]: column with bad row index
1.374 + stats [5]: bad row index
1.375 + stats [6]: n_row, # of rows of matrx
1.376 +
1.377 + -10 (unused; see symamd.c)
1.378 +
1.379 + -999 (unused; see symamd.c)
1.380 +
1.381 + Future versions may return more statistics in the stats array.
1.382 +
1.383 + Example:
1.384 +
1.385 + See http://www.cise.ufl.edu/research/sparse/colamd/example.c
1.386 + for a complete example.
1.387 +
1.388 + To order the columns of a 5-by-4 matrix with 11 nonzero entries in
1.389 + the following nonzero pattern
1.390 +
1.391 + x 0 x 0
1.392 + x 0 x x
1.393 + 0 x x 0
1.394 + 0 0 x x
1.395 + x x 0 0
1.396 +
1.397 + with default knobs and no output statistics, do the following:
1.398 +
1.399 + #include "colamd.h"
1.400 + #define ALEN 100
1.401 + int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
1.402 + int p [ ] = {0, 3, 5, 9, 11} ;
1.403 + int stats [COLAMD_STATS] ;
1.404 + colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
1.405 +
1.406 + The permutation is returned in the array p, and A is destroyed.
1.407 +
1.408 + ----------------------------------------------------------------------------
1.409 + symamd:
1.410 + ----------------------------------------------------------------------------
1.411 +
1.412 + C syntax:
1.413 +
1.414 + #include "colamd.h"
1.415 + int symamd (int n, int *A, int *p, int *perm,
1.416 + double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
1.417 + void (*allocate) (size_t, size_t), void (*release) (void *)) ;
1.418 + UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
1.419 + double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
1.420 + void (*allocate) (size_t, size_t), void (*release) (void *)) ;
1.421 +
1.422 + Purpose:
1.423 +
1.424 + The symamd routine computes an ordering P of a symmetric sparse
1.425 + matrix A such that the Cholesky factorization PAP' = LL' remains
1.426 + sparse. It is based on a column ordering of a matrix M constructed
1.427 + so that the nonzero pattern of M'M is the same as A. The matrix A
1.428 + is assumed to be symmetric; only the strictly lower triangular part
1.429 + is accessed. You must pass your selected memory allocator (usually
1.430 + calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
1.431 + memory for the temporary matrix M.
1.432 +
1.433 + Returns:
1.434 +
1.435 + TRUE (1) if successful, FALSE (0) otherwise.
1.436 +
1.437 + Arguments:
1.438 +
1.439 + int n ; Input argument.
1.440 +
1.441 + Number of rows and columns in the symmetrix matrix A.
1.442 + Restriction: n >= 0.
1.443 + Symamd returns FALSE if n is negative.
1.444 +
1.445 + int A [nnz] ; Input argument.
1.446 +
1.447 + A is an integer array of size nnz, where nnz = p [n].
1.448 +
1.449 + The row indices of the entries in column c of the matrix are
1.450 + held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a
1.451 + given column c need not be in ascending order, and duplicate
1.452 + row indices may be present. However, symamd will run faster
1.453 + if the columns are in sorted order with no duplicate entries.
1.454 +
1.455 + The matrix is 0-based. That is, rows are in the range 0 to
1.456 + n-1, and columns are in the range 0 to n-1. Symamd
1.457 + returns FALSE if any row index is out of range.
1.458 +
1.459 + The contents of A are not modified.
1.460 +
1.461 + int p [n+1] ; Input argument.
1.462 +
1.463 + p is an integer array of size n+1. On input, it holds the
1.464 + "pointers" for the column form of the matrix A. Column c of
1.465 + the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
1.466 + entry, p [0], must be zero, and p [c] <= p [c+1] must hold
1.467 + for all c in the range 0 to n-1. The value p [n] is
1.468 + thus the total number of entries in the pattern of the matrix A.
1.469 + Symamd returns FALSE if these conditions are not met.
1.470 +
1.471 + The contents of p are not modified.
1.472 +
1.473 + int perm [n+1] ; Output argument.
1.474 +
1.475 + On output, if symamd returns TRUE, the array perm holds the
1.476 + permutation P, where perm [0] is the first index in the new
1.477 + ordering, and perm [n-1] is the last. That is, perm [k] = j
1.478 + means that row and column j of A is the kth column in PAP',
1.479 + where k is in the range 0 to n-1 (perm [0] = j means
1.480 + that row and column j of A are the first row and column in
1.481 + PAP'). The array is used as a workspace during the ordering,
1.482 + which is why it must be of length n+1, not just n.
1.483 +
1.484 + double knobs [COLAMD_KNOBS] ; Input argument.
1.485 +
1.486 + See colamd_set_defaults for a description.
1.487 +
1.488 + int stats [COLAMD_STATS] ; Output argument.
1.489 +
1.490 + Statistics on the ordering, and error status.
1.491 + See colamd.h for related definitions.
1.492 + Symamd returns FALSE if stats is not present.
1.493 +
1.494 + stats [0]: number of dense or empty row and columns ignored
1.495 + (and ordered last in the output permutation
1.496 + perm). Note that a row/column can become
1.497 + "empty" if it contains only "dense" and/or
1.498 + "empty" columns/rows.
1.499 +
1.500 + stats [1]: (same as stats [0])
1.501 +
1.502 + stats [2]: number of garbage collections performed.
1.503 +
1.504 + stats [3]: status code. < 0 is an error code.
1.505 + > 1 is a warning or notice.
1.506 +
1.507 + 0 OK. Each column of the input matrix contained
1.508 + row indices in increasing order, with no
1.509 + duplicates.
1.510 +
1.511 + 1 OK, but columns of input matrix were jumbled
1.512 + (unsorted columns or duplicate entries). Symamd
1.513 + had to do some extra work to sort the matrix
1.514 + first and remove duplicate entries, but it
1.515 + still was able to return a valid permutation
1.516 + (return value of symamd was TRUE).
1.517 +
1.518 + stats [4]: highest numbered column that
1.519 + is unsorted or has duplicate
1.520 + entries.
1.521 + stats [5]: last seen duplicate or
1.522 + unsorted row index.
1.523 + stats [6]: number of duplicate or
1.524 + unsorted row indices.
1.525 +
1.526 + -1 A is a null pointer
1.527 +
1.528 + -2 p is a null pointer
1.529 +
1.530 + -3 (unused, see colamd.c)
1.531 +
1.532 + -4 n is negative
1.533 +
1.534 + stats [4]: n
1.535 +
1.536 + -5 number of nonzeros in matrix is negative
1.537 +
1.538 + stats [4]: # of nonzeros (p [n]).
1.539 +
1.540 + -6 p [0] is nonzero
1.541 +
1.542 + stats [4]: p [0]
1.543 +
1.544 + -7 (unused)
1.545 +
1.546 + -8 a column has a negative number of entries
1.547 +
1.548 + stats [4]: column with < 0 entries
1.549 + stats [5]: number of entries in col
1.550 +
1.551 + -9 a row index is out of bounds
1.552 +
1.553 + stats [4]: column with bad row index
1.554 + stats [5]: bad row index
1.555 + stats [6]: n_row, # of rows of matrx
1.556 +
1.557 + -10 out of memory (unable to allocate temporary
1.558 + workspace for M or count arrays using the
1.559 + "allocate" routine passed into symamd).
1.560 +
1.561 + Future versions may return more statistics in the stats array.
1.562 +
1.563 + void * (*allocate) (size_t, size_t)
1.564 +
1.565 + A pointer to a function providing memory allocation. The
1.566 + allocated memory must be returned initialized to zero. For a
1.567 + C application, this argument should normally be a pointer to
1.568 + calloc. For a MATLAB mexFunction, the routine mxCalloc is
1.569 + passed instead.
1.570 +
1.571 + void (*release) (size_t, size_t)
1.572 +
1.573 + A pointer to a function that frees memory allocated by the
1.574 + memory allocation routine above. For a C application, this
1.575 + argument should normally be a pointer to free. For a MATLAB
1.576 + mexFunction, the routine mxFree is passed instead.
1.577 +
1.578 +
1.579 + ----------------------------------------------------------------------------
1.580 + colamd_report:
1.581 + ----------------------------------------------------------------------------
1.582 +
1.583 + C syntax:
1.584 +
1.585 + #include "colamd.h"
1.586 + colamd_report (int stats [COLAMD_STATS]) ;
1.587 + colamd_l_report (UF_long stats [COLAMD_STATS]) ;
1.588 +
1.589 + Purpose:
1.590 +
1.591 + Prints the error status and statistics recorded in the stats
1.592 + array on the standard error output (for a standard C routine)
1.593 + or on the MATLAB output (for a mexFunction).
1.594 +
1.595 + Arguments:
1.596 +
1.597 + int stats [COLAMD_STATS] ; Input only. Statistics from colamd.
1.598 +
1.599 +
1.600 + ----------------------------------------------------------------------------
1.601 + symamd_report:
1.602 + ----------------------------------------------------------------------------
1.603 +
1.604 + C syntax:
1.605 +
1.606 + #include "colamd.h"
1.607 + symamd_report (int stats [COLAMD_STATS]) ;
1.608 + symamd_l_report (UF_long stats [COLAMD_STATS]) ;
1.609 +
1.610 + Purpose:
1.611 +
1.612 + Prints the error status and statistics recorded in the stats
1.613 + array on the standard error output (for a standard C routine)
1.614 + or on the MATLAB output (for a mexFunction).
1.615 +
1.616 + Arguments:
1.617 +
1.618 + int stats [COLAMD_STATS] ; Input only. Statistics from symamd.
1.619 +
1.620 +
1.621 +*/
1.622 +
1.623 +/* ========================================================================== */
1.624 +/* === Scaffolding code definitions ======================================== */
1.625 +/* ========================================================================== */
1.626 +
1.627 +/* Ensure that debugging is turned off: */
1.628 +#ifndef NDEBUG
1.629 +#define NDEBUG
1.630 +#endif
1.631 +
1.632 +/* turn on debugging by uncommenting the following line
1.633 + #undef NDEBUG
1.634 +*/
1.635 +
1.636 +/*
1.637 + Our "scaffolding code" philosophy: In our opinion, well-written library
1.638 + code should keep its "debugging" code, and just normally have it turned off
1.639 + by the compiler so as not to interfere with performance. This serves
1.640 + several purposes:
1.641 +
1.642 + (1) assertions act as comments to the reader, telling you what the code
1.643 + expects at that point. All assertions will always be true (unless
1.644 + there really is a bug, of course).
1.645 +
1.646 + (2) leaving in the scaffolding code assists anyone who would like to modify
1.647 + the code, or understand the algorithm (by reading the debugging output,
1.648 + one can get a glimpse into what the code is doing).
1.649 +
1.650 + (3) (gasp!) for actually finding bugs. This code has been heavily tested
1.651 + and "should" be fully functional and bug-free ... but you never know...
1.652 +
1.653 + The code will become outrageously slow when debugging is
1.654 + enabled. To control the level of debugging output, set an environment
1.655 + variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging,
1.656 + you should see the following message on the standard output:
1.657 +
1.658 + colamd: debug version, D = 1 (THIS WILL BE SLOW!)
1.659 +
1.660 + or a similar message for symamd. If you don't, then debugging has not
1.661 + been enabled.
1.662 +
1.663 +*/
1.664 +
1.665 +/* ========================================================================== */
1.666 +/* === Include files ======================================================== */
1.667 +/* ========================================================================== */
1.668 +
1.669 +#include "colamd.h"
1.670 +
1.671 +#if 0 /* by mao */
1.672 +#include <limits.h>
1.673 +#include <math.h>
1.674 +
1.675 +#ifdef MATLAB_MEX_FILE
1.676 +#include "mex.h"
1.677 +#include "matrix.h"
1.678 +#endif /* MATLAB_MEX_FILE */
1.679 +
1.680 +#if !defined (NPRINT) || !defined (NDEBUG)
1.681 +#include <stdio.h>
1.682 +#endif
1.683 +
1.684 +#ifndef NULL
1.685 +#define NULL ((void *) 0)
1.686 +#endif
1.687 +#endif
1.688 +
1.689 +/* ========================================================================== */
1.690 +/* === int or UF_long ======================================================= */
1.691 +/* ========================================================================== */
1.692 +
1.693 +#if 0 /* by mao */
1.694 +/* define UF_long */
1.695 +#include "UFconfig.h"
1.696 +#endif
1.697 +
1.698 +#ifdef DLONG
1.699 +
1.700 +#define Int UF_long
1.701 +#define ID UF_long_id
1.702 +#define Int_MAX UF_long_max
1.703 +
1.704 +#define COLAMD_recommended colamd_l_recommended
1.705 +#define COLAMD_set_defaults colamd_l_set_defaults
1.706 +#define COLAMD_MAIN colamd_l
1.707 +#define SYMAMD_MAIN symamd_l
1.708 +#define COLAMD_report colamd_l_report
1.709 +#define SYMAMD_report symamd_l_report
1.710 +
1.711 +#else
1.712 +
1.713 +#define Int int
1.714 +#define ID "%d"
1.715 +#define Int_MAX INT_MAX
1.716 +
1.717 +#define COLAMD_recommended colamd_recommended
1.718 +#define COLAMD_set_defaults colamd_set_defaults
1.719 +#define COLAMD_MAIN colamd
1.720 +#define SYMAMD_MAIN symamd
1.721 +#define COLAMD_report colamd_report
1.722 +#define SYMAMD_report symamd_report
1.723 +
1.724 +#endif
1.725 +
1.726 +/* ========================================================================== */
1.727 +/* === Row and Column structures ============================================ */
1.728 +/* ========================================================================== */
1.729 +
1.730 +/* User code that makes use of the colamd/symamd routines need not directly */
1.731 +/* reference these structures. They are used only for colamd_recommended. */
1.732 +
1.733 +typedef struct Colamd_Col_struct
1.734 +{
1.735 + Int start ; /* index for A of first row in this column, or DEAD */
1.736 + /* if column is dead */
1.737 + Int length ; /* number of rows in this column */
1.738 + union
1.739 + {
1.740 + Int thickness ; /* number of original columns represented by this */
1.741 + /* col, if the column is alive */
1.742 + Int parent ; /* parent in parent tree super-column structure, if */
1.743 + /* the column is dead */
1.744 + } shared1 ;
1.745 + union
1.746 + {
1.747 + Int score ; /* the score used to maintain heap, if col is alive */
1.748 + Int order ; /* pivot ordering of this column, if col is dead */
1.749 + } shared2 ;
1.750 + union
1.751 + {
1.752 + Int headhash ; /* head of a hash bucket, if col is at the head of */
1.753 + /* a degree list */
1.754 + Int hash ; /* hash value, if col is not in a degree list */
1.755 + Int prev ; /* previous column in degree list, if col is in a */
1.756 + /* degree list (but not at the head of a degree list) */
1.757 + } shared3 ;
1.758 + union
1.759 + {
1.760 + Int degree_next ; /* next column, if col is in a degree list */
1.761 + Int hash_next ; /* next column, if col is in a hash list */
1.762 + } shared4 ;
1.763 +
1.764 +} Colamd_Col ;
1.765 +
1.766 +typedef struct Colamd_Row_struct
1.767 +{
1.768 + Int start ; /* index for A of first col in this row */
1.769 + Int length ; /* number of principal columns in this row */
1.770 + union
1.771 + {
1.772 + Int degree ; /* number of principal & non-principal columns in row */
1.773 + Int p ; /* used as a row pointer in init_rows_cols () */
1.774 + } shared1 ;
1.775 + union
1.776 + {
1.777 + Int mark ; /* for computing set differences and marking dead rows*/
1.778 + Int first_column ;/* first column in row (used in garbage collection) */
1.779 + } shared2 ;
1.780 +
1.781 +} Colamd_Row ;
1.782 +
1.783 +/* ========================================================================== */
1.784 +/* === Definitions ========================================================== */
1.785 +/* ========================================================================== */
1.786 +
1.787 +/* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
1.788 +#define PUBLIC
1.789 +#define PRIVATE static
1.790 +
1.791 +#define DENSE_DEGREE(alpha,n) \
1.792 + ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
1.793 +
1.794 +#define MAX(a,b) (((a) > (b)) ? (a) : (b))
1.795 +#define MIN(a,b) (((a) < (b)) ? (a) : (b))
1.796 +
1.797 +#define ONES_COMPLEMENT(r) (-(r)-1)
1.798 +
1.799 +/* -------------------------------------------------------------------------- */
1.800 +/* Change for version 2.1: define TRUE and FALSE only if not yet defined */
1.801 +/* -------------------------------------------------------------------------- */
1.802 +
1.803 +#ifndef TRUE
1.804 +#define TRUE (1)
1.805 +#endif
1.806 +
1.807 +#ifndef FALSE
1.808 +#define FALSE (0)
1.809 +#endif
1.810 +
1.811 +/* -------------------------------------------------------------------------- */
1.812 +
1.813 +#define EMPTY (-1)
1.814 +
1.815 +/* Row and column status */
1.816 +#define ALIVE (0)
1.817 +#define DEAD (-1)
1.818 +
1.819 +/* Column status */
1.820 +#define DEAD_PRINCIPAL (-1)
1.821 +#define DEAD_NON_PRINCIPAL (-2)
1.822 +
1.823 +/* Macros for row and column status update and checking. */
1.824 +#define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
1.825 +#define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
1.826 +#define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
1.827 +#define COL_IS_DEAD(c) (Col [c].start < ALIVE)
1.828 +#define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
1.829 +#define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
1.830 +#define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
1.831 +#define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
1.832 +#define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
1.833 +
1.834 +/* ========================================================================== */
1.835 +/* === Colamd reporting mechanism =========================================== */
1.836 +/* ========================================================================== */
1.837 +
1.838 +#if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
1.839 +/* In MATLAB, matrices are 1-based to the user, but 0-based internally */
1.840 +#define INDEX(i) ((i)+1)
1.841 +#else
1.842 +/* In C, matrices are 0-based and indices are reported as such in *_report */
1.843 +#define INDEX(i) (i)
1.844 +#endif
1.845 +
1.846 +/* All output goes through the PRINTF macro. */
1.847 +#define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
1.848 +
1.849 +/* ========================================================================== */
1.850 +/* === Prototypes of PRIVATE routines ======================================= */
1.851 +/* ========================================================================== */
1.852 +
1.853 +PRIVATE Int init_rows_cols
1.854 +(
1.855 + Int n_row,
1.856 + Int n_col,
1.857 + Colamd_Row Row [],
1.858 + Colamd_Col Col [],
1.859 + Int A [],
1.860 + Int p [],
1.861 + Int stats [COLAMD_STATS]
1.862 +) ;
1.863 +
1.864 +PRIVATE void init_scoring
1.865 +(
1.866 + Int n_row,
1.867 + Int n_col,
1.868 + Colamd_Row Row [],
1.869 + Colamd_Col Col [],
1.870 + Int A [],
1.871 + Int head [],
1.872 + double knobs [COLAMD_KNOBS],
1.873 + Int *p_n_row2,
1.874 + Int *p_n_col2,
1.875 + Int *p_max_deg
1.876 +) ;
1.877 +
1.878 +PRIVATE Int find_ordering
1.879 +(
1.880 + Int n_row,
1.881 + Int n_col,
1.882 + Int Alen,
1.883 + Colamd_Row Row [],
1.884 + Colamd_Col Col [],
1.885 + Int A [],
1.886 + Int head [],
1.887 + Int n_col2,
1.888 + Int max_deg,
1.889 + Int pfree,
1.890 + Int aggressive
1.891 +) ;
1.892 +
1.893 +PRIVATE void order_children
1.894 +(
1.895 + Int n_col,
1.896 + Colamd_Col Col [],
1.897 + Int p []
1.898 +) ;
1.899 +
1.900 +PRIVATE void detect_super_cols
1.901 +(
1.902 +
1.903 +#ifndef NDEBUG
1.904 + Int n_col,
1.905 + Colamd_Row Row [],
1.906 +#endif /* NDEBUG */
1.907 +
1.908 + Colamd_Col Col [],
1.909 + Int A [],
1.910 + Int head [],
1.911 + Int row_start,
1.912 + Int row_length
1.913 +) ;
1.914 +
1.915 +PRIVATE Int garbage_collection
1.916 +(
1.917 + Int n_row,
1.918 + Int n_col,
1.919 + Colamd_Row Row [],
1.920 + Colamd_Col Col [],
1.921 + Int A [],
1.922 + Int *pfree
1.923 +) ;
1.924 +
1.925 +PRIVATE Int clear_mark
1.926 +(
1.927 + Int tag_mark,
1.928 + Int max_mark,
1.929 + Int n_row,
1.930 + Colamd_Row Row []
1.931 +) ;
1.932 +
1.933 +PRIVATE void print_report
1.934 +(
1.935 + char *method,
1.936 + Int stats [COLAMD_STATS]
1.937 +) ;
1.938 +
1.939 +/* ========================================================================== */
1.940 +/* === Debugging prototypes and definitions ================================= */
1.941 +/* ========================================================================== */
1.942 +
1.943 +#ifndef NDEBUG
1.944 +
1.945 +#if 0 /* by mao */
1.946 +#include <assert.h>
1.947 +#endif
1.948 +
1.949 +/* colamd_debug is the *ONLY* global variable, and is only */
1.950 +/* present when debugging */
1.951 +
1.952 +PRIVATE Int colamd_debug = 0 ; /* debug print level */
1.953 +
1.954 +#define DEBUG0(params) { PRINTF (params) ; }
1.955 +#define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
1.956 +#define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
1.957 +#define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
1.958 +#define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
1.959 +
1.960 +#if 0 /* by mao */
1.961 +#ifdef MATLAB_MEX_FILE
1.962 +#define ASSERT(expression) (mxAssert ((expression), ""))
1.963 +#else
1.964 +#define ASSERT(expression) (assert (expression))
1.965 +#endif /* MATLAB_MEX_FILE */
1.966 +#else
1.967 +#define ASSERT xassert
1.968 +#endif
1.969 +
1.970 +PRIVATE void colamd_get_debug /* gets the debug print level from getenv */
1.971 +(
1.972 + char *method
1.973 +) ;
1.974 +
1.975 +PRIVATE void debug_deg_lists
1.976 +(
1.977 + Int n_row,
1.978 + Int n_col,
1.979 + Colamd_Row Row [],
1.980 + Colamd_Col Col [],
1.981 + Int head [],
1.982 + Int min_score,
1.983 + Int should,
1.984 + Int max_deg
1.985 +) ;
1.986 +
1.987 +PRIVATE void debug_mark
1.988 +(
1.989 + Int n_row,
1.990 + Colamd_Row Row [],
1.991 + Int tag_mark,
1.992 + Int max_mark
1.993 +) ;
1.994 +
1.995 +PRIVATE void debug_matrix
1.996 +(
1.997 + Int n_row,
1.998 + Int n_col,
1.999 + Colamd_Row Row [],
1.1000 + Colamd_Col Col [],
1.1001 + Int A []
1.1002 +) ;
1.1003 +
1.1004 +PRIVATE void debug_structures
1.1005 +(
1.1006 + Int n_row,
1.1007 + Int n_col,
1.1008 + Colamd_Row Row [],
1.1009 + Colamd_Col Col [],
1.1010 + Int A [],
1.1011 + Int n_col2
1.1012 +) ;
1.1013 +
1.1014 +#else /* NDEBUG */
1.1015 +
1.1016 +/* === No debugging ========================================================= */
1.1017 +
1.1018 +#define DEBUG0(params) ;
1.1019 +#define DEBUG1(params) ;
1.1020 +#define DEBUG2(params) ;
1.1021 +#define DEBUG3(params) ;
1.1022 +#define DEBUG4(params) ;
1.1023 +
1.1024 +#define ASSERT(expression)
1.1025 +
1.1026 +#endif /* NDEBUG */
1.1027 +
1.1028 +/* ========================================================================== */
1.1029 +/* === USER-CALLABLE ROUTINES: ============================================== */
1.1030 +/* ========================================================================== */
1.1031 +
1.1032 +/* ========================================================================== */
1.1033 +/* === colamd_recommended =================================================== */
1.1034 +/* ========================================================================== */
1.1035 +
1.1036 +/*
1.1037 + The colamd_recommended routine returns the suggested size for Alen. This
1.1038 + value has been determined to provide good balance between the number of
1.1039 + garbage collections and the memory requirements for colamd. If any
1.1040 + argument is negative, or if integer overflow occurs, a 0 is returned as an
1.1041 + error condition. 2*nnz space is required for the row and column
1.1042 + indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
1.1043 + required for the Col and Row arrays, respectively, which are internal to
1.1044 + colamd (roughly 6*n_col + 4*n_row). An additional n_col space is the
1.1045 + minimal amount of "elbow room", and nnz/5 more space is recommended for
1.1046 + run time efficiency.
1.1047 +
1.1048 + Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
1.1049 +
1.1050 + This function is not needed when using symamd.
1.1051 +*/
1.1052 +
1.1053 +/* add two values of type size_t, and check for integer overflow */
1.1054 +static size_t t_add (size_t a, size_t b, int *ok)
1.1055 +{
1.1056 + (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
1.1057 + return ((*ok) ? (a + b) : 0) ;
1.1058 +}
1.1059 +
1.1060 +/* compute a*k where k is a small integer, and check for integer overflow */
1.1061 +static size_t t_mult (size_t a, size_t k, int *ok)
1.1062 +{
1.1063 + size_t i, s = 0 ;
1.1064 + for (i = 0 ; i < k ; i++)
1.1065 + {
1.1066 + s = t_add (s, a, ok) ;
1.1067 + }
1.1068 + return (s) ;
1.1069 +}
1.1070 +
1.1071 +/* size of the Col and Row structures */
1.1072 +#define COLAMD_C(n_col,ok) \
1.1073 + ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
1.1074 +
1.1075 +#define COLAMD_R(n_row,ok) \
1.1076 + ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
1.1077 +
1.1078 +
1.1079 +PUBLIC size_t COLAMD_recommended /* returns recommended value of Alen. */
1.1080 +(
1.1081 + /* === Parameters ======================================================= */
1.1082 +
1.1083 + Int nnz, /* number of nonzeros in A */
1.1084 + Int n_row, /* number of rows in A */
1.1085 + Int n_col /* number of columns in A */
1.1086 +)
1.1087 +{
1.1088 + size_t s, c, r ;
1.1089 + int ok = TRUE ;
1.1090 + if (nnz < 0 || n_row < 0 || n_col < 0)
1.1091 + {
1.1092 + return (0) ;
1.1093 + }
1.1094 + s = t_mult (nnz, 2, &ok) ; /* 2*nnz */
1.1095 + c = COLAMD_C (n_col, &ok) ; /* size of column structures */
1.1096 + r = COLAMD_R (n_row, &ok) ; /* size of row structures */
1.1097 + s = t_add (s, c, &ok) ;
1.1098 + s = t_add (s, r, &ok) ;
1.1099 + s = t_add (s, n_col, &ok) ; /* elbow room */
1.1100 + s = t_add (s, nnz/5, &ok) ; /* elbow room */
1.1101 + ok = ok && (s < Int_MAX) ;
1.1102 + return (ok ? s : 0) ;
1.1103 +}
1.1104 +
1.1105 +
1.1106 +/* ========================================================================== */
1.1107 +/* === colamd_set_defaults ================================================== */
1.1108 +/* ========================================================================== */
1.1109 +
1.1110 +/*
1.1111 + The colamd_set_defaults routine sets the default values of the user-
1.1112 + controllable parameters for colamd and symamd:
1.1113 +
1.1114 + Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
1.1115 + entries are removed prior to ordering. Columns with more than
1.1116 + max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
1.1117 + prior to ordering, and placed last in the output column ordering.
1.1118 +
1.1119 + Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
1.1120 + entries are removed prior to ordering, and placed last in the
1.1121 + output ordering.
1.1122 +
1.1123 + knobs [0] dense row control
1.1124 +
1.1125 + knobs [1] dense column control
1.1126 +
1.1127 + knobs [2] if nonzero, do aggresive absorption
1.1128 +
1.1129 + knobs [3..19] unused, but future versions might use this
1.1130 +
1.1131 +*/
1.1132 +
1.1133 +PUBLIC void COLAMD_set_defaults
1.1134 +(
1.1135 + /* === Parameters ======================================================= */
1.1136 +
1.1137 + double knobs [COLAMD_KNOBS] /* knob array */
1.1138 +)
1.1139 +{
1.1140 + /* === Local variables ================================================== */
1.1141 +
1.1142 + Int i ;
1.1143 +
1.1144 + if (!knobs)
1.1145 + {
1.1146 + return ; /* no knobs to initialize */
1.1147 + }
1.1148 + for (i = 0 ; i < COLAMD_KNOBS ; i++)
1.1149 + {
1.1150 + knobs [i] = 0 ;
1.1151 + }
1.1152 + knobs [COLAMD_DENSE_ROW] = 10 ;
1.1153 + knobs [COLAMD_DENSE_COL] = 10 ;
1.1154 + knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
1.1155 +}
1.1156 +
1.1157 +
1.1158 +/* ========================================================================== */
1.1159 +/* === symamd =============================================================== */
1.1160 +/* ========================================================================== */
1.1161 +
1.1162 +PUBLIC Int SYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */
1.1163 +(
1.1164 + /* === Parameters ======================================================= */
1.1165 +
1.1166 + Int n, /* number of rows and columns of A */
1.1167 + Int A [], /* row indices of A */
1.1168 + Int p [], /* column pointers of A */
1.1169 + Int perm [], /* output permutation, size n+1 */
1.1170 + double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
1.1171 + Int stats [COLAMD_STATS], /* output statistics and error codes */
1.1172 + void * (*allocate) (size_t, size_t),
1.1173 + /* pointer to calloc (ANSI C) or */
1.1174 + /* mxCalloc (for MATLAB mexFunction) */
1.1175 + void (*release) (void *)
1.1176 + /* pointer to free (ANSI C) or */
1.1177 + /* mxFree (for MATLAB mexFunction) */
1.1178 +)
1.1179 +{
1.1180 + /* === Local variables ================================================== */
1.1181 +
1.1182 + Int *count ; /* length of each column of M, and col pointer*/
1.1183 + Int *mark ; /* mark array for finding duplicate entries */
1.1184 + Int *M ; /* row indices of matrix M */
1.1185 + size_t Mlen ; /* length of M */
1.1186 + Int n_row ; /* number of rows in M */
1.1187 + Int nnz ; /* number of entries in A */
1.1188 + Int i ; /* row index of A */
1.1189 + Int j ; /* column index of A */
1.1190 + Int k ; /* row index of M */
1.1191 + Int mnz ; /* number of nonzeros in M */
1.1192 + Int pp ; /* index into a column of A */
1.1193 + Int last_row ; /* last row seen in the current column */
1.1194 + Int length ; /* number of nonzeros in a column */
1.1195 +
1.1196 + double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */
1.1197 + double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */
1.1198 +
1.1199 +#ifndef NDEBUG
1.1200 + colamd_get_debug ("symamd") ;
1.1201 +#endif /* NDEBUG */
1.1202 +
1.1203 + /* === Check the input arguments ======================================== */
1.1204 +
1.1205 + if (!stats)
1.1206 + {
1.1207 + DEBUG0 (("symamd: stats not present\n")) ;
1.1208 + return (FALSE) ;
1.1209 + }
1.1210 + for (i = 0 ; i < COLAMD_STATS ; i++)
1.1211 + {
1.1212 + stats [i] = 0 ;
1.1213 + }
1.1214 + stats [COLAMD_STATUS] = COLAMD_OK ;
1.1215 + stats [COLAMD_INFO1] = -1 ;
1.1216 + stats [COLAMD_INFO2] = -1 ;
1.1217 +
1.1218 + if (!A)
1.1219 + {
1.1220 + stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1.1221 + DEBUG0 (("symamd: A not present\n")) ;
1.1222 + return (FALSE) ;
1.1223 + }
1.1224 +
1.1225 + if (!p) /* p is not present */
1.1226 + {
1.1227 + stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1.1228 + DEBUG0 (("symamd: p not present\n")) ;
1.1229 + return (FALSE) ;
1.1230 + }
1.1231 +
1.1232 + if (n < 0) /* n must be >= 0 */
1.1233 + {
1.1234 + stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1.1235 + stats [COLAMD_INFO1] = n ;
1.1236 + DEBUG0 (("symamd: n negative %d\n", n)) ;
1.1237 + return (FALSE) ;
1.1238 + }
1.1239 +
1.1240 + nnz = p [n] ;
1.1241 + if (nnz < 0) /* nnz must be >= 0 */
1.1242 + {
1.1243 + stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1.1244 + stats [COLAMD_INFO1] = nnz ;
1.1245 + DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
1.1246 + return (FALSE) ;
1.1247 + }
1.1248 +
1.1249 + if (p [0] != 0)
1.1250 + {
1.1251 + stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1.1252 + stats [COLAMD_INFO1] = p [0] ;
1.1253 + DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
1.1254 + return (FALSE) ;
1.1255 + }
1.1256 +
1.1257 + /* === If no knobs, set default knobs =================================== */
1.1258 +
1.1259 + if (!knobs)
1.1260 + {
1.1261 + COLAMD_set_defaults (default_knobs) ;
1.1262 + knobs = default_knobs ;
1.1263 + }
1.1264 +
1.1265 + /* === Allocate count and mark ========================================== */
1.1266 +
1.1267 + count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1.1268 + if (!count)
1.1269 + {
1.1270 + stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1.1271 + DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
1.1272 + return (FALSE) ;
1.1273 + }
1.1274 +
1.1275 + mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1.1276 + if (!mark)
1.1277 + {
1.1278 + stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1.1279 + (*release) ((void *) count) ;
1.1280 + DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
1.1281 + return (FALSE) ;
1.1282 + }
1.1283 +
1.1284 + /* === Compute column counts of M, check if A is valid ================== */
1.1285 +
1.1286 + stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
1.1287 +
1.1288 + for (i = 0 ; i < n ; i++)
1.1289 + {
1.1290 + mark [i] = -1 ;
1.1291 + }
1.1292 +
1.1293 + for (j = 0 ; j < n ; j++)
1.1294 + {
1.1295 + last_row = -1 ;
1.1296 +
1.1297 + length = p [j+1] - p [j] ;
1.1298 + if (length < 0)
1.1299 + {
1.1300 + /* column pointers must be non-decreasing */
1.1301 + stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1.1302 + stats [COLAMD_INFO1] = j ;
1.1303 + stats [COLAMD_INFO2] = length ;
1.1304 + (*release) ((void *) count) ;
1.1305 + (*release) ((void *) mark) ;
1.1306 + DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
1.1307 + return (FALSE) ;
1.1308 + }
1.1309 +
1.1310 + for (pp = p [j] ; pp < p [j+1] ; pp++)
1.1311 + {
1.1312 + i = A [pp] ;
1.1313 + if (i < 0 || i >= n)
1.1314 + {
1.1315 + /* row index i, in column j, is out of bounds */
1.1316 + stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1.1317 + stats [COLAMD_INFO1] = j ;
1.1318 + stats [COLAMD_INFO2] = i ;
1.1319 + stats [COLAMD_INFO3] = n ;
1.1320 + (*release) ((void *) count) ;
1.1321 + (*release) ((void *) mark) ;
1.1322 + DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
1.1323 + return (FALSE) ;
1.1324 + }
1.1325 +
1.1326 + if (i <= last_row || mark [i] == j)
1.1327 + {
1.1328 + /* row index is unsorted or repeated (or both), thus col */
1.1329 + /* is jumbled. This is a notice, not an error condition. */
1.1330 + stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1.1331 + stats [COLAMD_INFO1] = j ;
1.1332 + stats [COLAMD_INFO2] = i ;
1.1333 + (stats [COLAMD_INFO3]) ++ ;
1.1334 + DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
1.1335 + }
1.1336 +
1.1337 + if (i > j && mark [i] != j)
1.1338 + {
1.1339 + /* row k of M will contain column indices i and j */
1.1340 + count [i]++ ;
1.1341 + count [j]++ ;
1.1342 + }
1.1343 +
1.1344 + /* mark the row as having been seen in this column */
1.1345 + mark [i] = j ;
1.1346 +
1.1347 + last_row = i ;
1.1348 + }
1.1349 + }
1.1350 +
1.1351 + /* v2.4: removed free(mark) */
1.1352 +
1.1353 + /* === Compute column pointers of M ===================================== */
1.1354 +
1.1355 + /* use output permutation, perm, for column pointers of M */
1.1356 + perm [0] = 0 ;
1.1357 + for (j = 1 ; j <= n ; j++)
1.1358 + {
1.1359 + perm [j] = perm [j-1] + count [j-1] ;
1.1360 + }
1.1361 + for (j = 0 ; j < n ; j++)
1.1362 + {
1.1363 + count [j] = perm [j] ;
1.1364 + }
1.1365 +
1.1366 + /* === Construct M ====================================================== */
1.1367 +
1.1368 + mnz = perm [n] ;
1.1369 + n_row = mnz / 2 ;
1.1370 + Mlen = COLAMD_recommended (mnz, n_row, n) ;
1.1371 + M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
1.1372 + DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
1.1373 + n_row, n, mnz, (double) Mlen)) ;
1.1374 +
1.1375 + if (!M)
1.1376 + {
1.1377 + stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1.1378 + (*release) ((void *) count) ;
1.1379 + (*release) ((void *) mark) ;
1.1380 + DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
1.1381 + return (FALSE) ;
1.1382 + }
1.1383 +
1.1384 + k = 0 ;
1.1385 +
1.1386 + if (stats [COLAMD_STATUS] == COLAMD_OK)
1.1387 + {
1.1388 + /* Matrix is OK */
1.1389 + for (j = 0 ; j < n ; j++)
1.1390 + {
1.1391 + ASSERT (p [j+1] - p [j] >= 0) ;
1.1392 + for (pp = p [j] ; pp < p [j+1] ; pp++)
1.1393 + {
1.1394 + i = A [pp] ;
1.1395 + ASSERT (i >= 0 && i < n) ;
1.1396 + if (i > j)
1.1397 + {
1.1398 + /* row k of M contains column indices i and j */
1.1399 + M [count [i]++] = k ;
1.1400 + M [count [j]++] = k ;
1.1401 + k++ ;
1.1402 + }
1.1403 + }
1.1404 + }
1.1405 + }
1.1406 + else
1.1407 + {
1.1408 + /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */
1.1409 + DEBUG0 (("symamd: Duplicates in A.\n")) ;
1.1410 + for (i = 0 ; i < n ; i++)
1.1411 + {
1.1412 + mark [i] = -1 ;
1.1413 + }
1.1414 + for (j = 0 ; j < n ; j++)
1.1415 + {
1.1416 + ASSERT (p [j+1] - p [j] >= 0) ;
1.1417 + for (pp = p [j] ; pp < p [j+1] ; pp++)
1.1418 + {
1.1419 + i = A [pp] ;
1.1420 + ASSERT (i >= 0 && i < n) ;
1.1421 + if (i > j && mark [i] != j)
1.1422 + {
1.1423 + /* row k of M contains column indices i and j */
1.1424 + M [count [i]++] = k ;
1.1425 + M [count [j]++] = k ;
1.1426 + k++ ;
1.1427 + mark [i] = j ;
1.1428 + }
1.1429 + }
1.1430 + }
1.1431 + /* v2.4: free(mark) moved below */
1.1432 + }
1.1433 +
1.1434 + /* count and mark no longer needed */
1.1435 + (*release) ((void *) count) ;
1.1436 + (*release) ((void *) mark) ; /* v2.4: free (mark) moved here */
1.1437 + ASSERT (k == n_row) ;
1.1438 +
1.1439 + /* === Adjust the knobs for M =========================================== */
1.1440 +
1.1441 + for (i = 0 ; i < COLAMD_KNOBS ; i++)
1.1442 + {
1.1443 + cknobs [i] = knobs [i] ;
1.1444 + }
1.1445 +
1.1446 + /* there are no dense rows in M */
1.1447 + cknobs [COLAMD_DENSE_ROW] = -1 ;
1.1448 + cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
1.1449 +
1.1450 + /* === Order the columns of M =========================================== */
1.1451 +
1.1452 + /* v2.4: colamd cannot fail here, so the error check is removed */
1.1453 + (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
1.1454 +
1.1455 + /* Note that the output permutation is now in perm */
1.1456 +
1.1457 + /* === get the statistics for symamd from colamd ======================== */
1.1458 +
1.1459 + /* a dense column in colamd means a dense row and col in symamd */
1.1460 + stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
1.1461 +
1.1462 + /* === Free M =========================================================== */
1.1463 +
1.1464 + (*release) ((void *) M) ;
1.1465 + DEBUG0 (("symamd: done.\n")) ;
1.1466 + return (TRUE) ;
1.1467 +
1.1468 +}
1.1469 +
1.1470 +/* ========================================================================== */
1.1471 +/* === colamd =============================================================== */
1.1472 +/* ========================================================================== */
1.1473 +
1.1474 +/*
1.1475 + The colamd routine computes a column ordering Q of a sparse matrix
1.1476 + A such that the LU factorization P(AQ) = LU remains sparse, where P is
1.1477 + selected via partial pivoting. The routine can also be viewed as
1.1478 + providing a permutation Q such that the Cholesky factorization
1.1479 + (AQ)'(AQ) = LL' remains sparse.
1.1480 +*/
1.1481 +
1.1482 +PUBLIC Int COLAMD_MAIN /* returns TRUE if successful, FALSE otherwise*/
1.1483 +(
1.1484 + /* === Parameters ======================================================= */
1.1485 +
1.1486 + Int n_row, /* number of rows in A */
1.1487 + Int n_col, /* number of columns in A */
1.1488 + Int Alen, /* length of A */
1.1489 + Int A [], /* row indices of A */
1.1490 + Int p [], /* pointers to columns in A */
1.1491 + double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
1.1492 + Int stats [COLAMD_STATS] /* output statistics and error codes */
1.1493 +)
1.1494 +{
1.1495 + /* === Local variables ================================================== */
1.1496 +
1.1497 + Int i ; /* loop index */
1.1498 + Int nnz ; /* nonzeros in A */
1.1499 + size_t Row_size ; /* size of Row [], in integers */
1.1500 + size_t Col_size ; /* size of Col [], in integers */
1.1501 + size_t need ; /* minimum required length of A */
1.1502 + Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */
1.1503 + Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */
1.1504 + Int n_col2 ; /* number of non-dense, non-empty columns */
1.1505 + Int n_row2 ; /* number of non-dense, non-empty rows */
1.1506 + Int ngarbage ; /* number of garbage collections performed */
1.1507 + Int max_deg ; /* maximum row degree */
1.1508 + double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
1.1509 + Int aggressive ; /* do aggressive absorption */
1.1510 + int ok ;
1.1511 +
1.1512 +#ifndef NDEBUG
1.1513 + colamd_get_debug ("colamd") ;
1.1514 +#endif /* NDEBUG */
1.1515 +
1.1516 + /* === Check the input arguments ======================================== */
1.1517 +
1.1518 + if (!stats)
1.1519 + {
1.1520 + DEBUG0 (("colamd: stats not present\n")) ;
1.1521 + return (FALSE) ;
1.1522 + }
1.1523 + for (i = 0 ; i < COLAMD_STATS ; i++)
1.1524 + {
1.1525 + stats [i] = 0 ;
1.1526 + }
1.1527 + stats [COLAMD_STATUS] = COLAMD_OK ;
1.1528 + stats [COLAMD_INFO1] = -1 ;
1.1529 + stats [COLAMD_INFO2] = -1 ;
1.1530 +
1.1531 + if (!A) /* A is not present */
1.1532 + {
1.1533 + stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1.1534 + DEBUG0 (("colamd: A not present\n")) ;
1.1535 + return (FALSE) ;
1.1536 + }
1.1537 +
1.1538 + if (!p) /* p is not present */
1.1539 + {
1.1540 + stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1.1541 + DEBUG0 (("colamd: p not present\n")) ;
1.1542 + return (FALSE) ;
1.1543 + }
1.1544 +
1.1545 + if (n_row < 0) /* n_row must be >= 0 */
1.1546 + {
1.1547 + stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
1.1548 + stats [COLAMD_INFO1] = n_row ;
1.1549 + DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
1.1550 + return (FALSE) ;
1.1551 + }
1.1552 +
1.1553 + if (n_col < 0) /* n_col must be >= 0 */
1.1554 + {
1.1555 + stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1.1556 + stats [COLAMD_INFO1] = n_col ;
1.1557 + DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
1.1558 + return (FALSE) ;
1.1559 + }
1.1560 +
1.1561 + nnz = p [n_col] ;
1.1562 + if (nnz < 0) /* nnz must be >= 0 */
1.1563 + {
1.1564 + stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1.1565 + stats [COLAMD_INFO1] = nnz ;
1.1566 + DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
1.1567 + return (FALSE) ;
1.1568 + }
1.1569 +
1.1570 + if (p [0] != 0)
1.1571 + {
1.1572 + stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1.1573 + stats [COLAMD_INFO1] = p [0] ;
1.1574 + DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
1.1575 + return (FALSE) ;
1.1576 + }
1.1577 +
1.1578 + /* === If no knobs, set default knobs =================================== */
1.1579 +
1.1580 + if (!knobs)
1.1581 + {
1.1582 + COLAMD_set_defaults (default_knobs) ;
1.1583 + knobs = default_knobs ;
1.1584 + }
1.1585 +
1.1586 + aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
1.1587 +
1.1588 + /* === Allocate the Row and Col arrays from array A ===================== */
1.1589 +
1.1590 + ok = TRUE ;
1.1591 + Col_size = COLAMD_C (n_col, &ok) ; /* size of Col array of structs */
1.1592 + Row_size = COLAMD_R (n_row, &ok) ; /* size of Row array of structs */
1.1593 +
1.1594 + /* need = 2*nnz + n_col + Col_size + Row_size ; */
1.1595 + need = t_mult (nnz, 2, &ok) ;
1.1596 + need = t_add (need, n_col, &ok) ;
1.1597 + need = t_add (need, Col_size, &ok) ;
1.1598 + need = t_add (need, Row_size, &ok) ;
1.1599 +
1.1600 + if (!ok || need > (size_t) Alen || need > Int_MAX)
1.1601 + {
1.1602 + /* not enough space in array A to perform the ordering */
1.1603 + stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
1.1604 + stats [COLAMD_INFO1] = need ;
1.1605 + stats [COLAMD_INFO2] = Alen ;
1.1606 + DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
1.1607 + return (FALSE) ;
1.1608 + }
1.1609 +
1.1610 + Alen -= Col_size + Row_size ;
1.1611 + Col = (Colamd_Col *) &A [Alen] ;
1.1612 + Row = (Colamd_Row *) &A [Alen + Col_size] ;
1.1613 +
1.1614 + /* === Construct the row and column data structures ===================== */
1.1615 +
1.1616 + if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
1.1617 + {
1.1618 + /* input matrix is invalid */
1.1619 + DEBUG0 (("colamd: Matrix invalid\n")) ;
1.1620 + return (FALSE) ;
1.1621 + }
1.1622 +
1.1623 + /* === Initialize scores, kill dense rows/columns ======================= */
1.1624 +
1.1625 + init_scoring (n_row, n_col, Row, Col, A, p, knobs,
1.1626 + &n_row2, &n_col2, &max_deg) ;
1.1627 +
1.1628 + /* === Order the supercolumns =========================================== */
1.1629 +
1.1630 + ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
1.1631 + n_col2, max_deg, 2*nnz, aggressive) ;
1.1632 +
1.1633 + /* === Order the non-principal columns ================================== */
1.1634 +
1.1635 + order_children (n_col, Col, p) ;
1.1636 +
1.1637 + /* === Return statistics in stats ======================================= */
1.1638 +
1.1639 + stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
1.1640 + stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
1.1641 + stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
1.1642 + DEBUG0 (("colamd: done.\n")) ;
1.1643 + return (TRUE) ;
1.1644 +}
1.1645 +
1.1646 +
1.1647 +/* ========================================================================== */
1.1648 +/* === colamd_report ======================================================== */
1.1649 +/* ========================================================================== */
1.1650 +
1.1651 +PUBLIC void COLAMD_report
1.1652 +(
1.1653 + Int stats [COLAMD_STATS]
1.1654 +)
1.1655 +{
1.1656 + print_report ("colamd", stats) ;
1.1657 +}
1.1658 +
1.1659 +
1.1660 +/* ========================================================================== */
1.1661 +/* === symamd_report ======================================================== */
1.1662 +/* ========================================================================== */
1.1663 +
1.1664 +PUBLIC void SYMAMD_report
1.1665 +(
1.1666 + Int stats [COLAMD_STATS]
1.1667 +)
1.1668 +{
1.1669 + print_report ("symamd", stats) ;
1.1670 +}
1.1671 +
1.1672 +
1.1673 +
1.1674 +/* ========================================================================== */
1.1675 +/* === NON-USER-CALLABLE ROUTINES: ========================================== */
1.1676 +/* ========================================================================== */
1.1677 +
1.1678 +/* There are no user-callable routines beyond this point in the file */
1.1679 +
1.1680 +
1.1681 +/* ========================================================================== */
1.1682 +/* === init_rows_cols ======================================================= */
1.1683 +/* ========================================================================== */
1.1684 +
1.1685 +/*
1.1686 + Takes the column form of the matrix in A and creates the row form of the
1.1687 + matrix. Also, row and column attributes are stored in the Col and Row
1.1688 + structs. If the columns are un-sorted or contain duplicate row indices,
1.1689 + this routine will also sort and remove duplicate row indices from the
1.1690 + column form of the matrix. Returns FALSE if the matrix is invalid,
1.1691 + TRUE otherwise. Not user-callable.
1.1692 +*/
1.1693 +
1.1694 +PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
1.1695 +(
1.1696 + /* === Parameters ======================================================= */
1.1697 +
1.1698 + Int n_row, /* number of rows of A */
1.1699 + Int n_col, /* number of columns of A */
1.1700 + Colamd_Row Row [], /* of size n_row+1 */
1.1701 + Colamd_Col Col [], /* of size n_col+1 */
1.1702 + Int A [], /* row indices of A, of size Alen */
1.1703 + Int p [], /* pointers to columns in A, of size n_col+1 */
1.1704 + Int stats [COLAMD_STATS] /* colamd statistics */
1.1705 +)
1.1706 +{
1.1707 + /* === Local variables ================================================== */
1.1708 +
1.1709 + Int col ; /* a column index */
1.1710 + Int row ; /* a row index */
1.1711 + Int *cp ; /* a column pointer */
1.1712 + Int *cp_end ; /* a pointer to the end of a column */
1.1713 + Int *rp ; /* a row pointer */
1.1714 + Int *rp_end ; /* a pointer to the end of a row */
1.1715 + Int last_row ; /* previous row */
1.1716 +
1.1717 + /* === Initialize columns, and check column pointers ==================== */
1.1718 +
1.1719 + for (col = 0 ; col < n_col ; col++)
1.1720 + {
1.1721 + Col [col].start = p [col] ;
1.1722 + Col [col].length = p [col+1] - p [col] ;
1.1723 +
1.1724 + if (Col [col].length < 0)
1.1725 + {
1.1726 + /* column pointers must be non-decreasing */
1.1727 + stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1.1728 + stats [COLAMD_INFO1] = col ;
1.1729 + stats [COLAMD_INFO2] = Col [col].length ;
1.1730 + DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
1.1731 + return (FALSE) ;
1.1732 + }
1.1733 +
1.1734 + Col [col].shared1.thickness = 1 ;
1.1735 + Col [col].shared2.score = 0 ;
1.1736 + Col [col].shared3.prev = EMPTY ;
1.1737 + Col [col].shared4.degree_next = EMPTY ;
1.1738 + }
1.1739 +
1.1740 + /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
1.1741 +
1.1742 + /* === Scan columns, compute row degrees, and check row indices ========= */
1.1743 +
1.1744 + stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
1.1745 +
1.1746 + for (row = 0 ; row < n_row ; row++)
1.1747 + {
1.1748 + Row [row].length = 0 ;
1.1749 + Row [row].shared2.mark = -1 ;
1.1750 + }
1.1751 +
1.1752 + for (col = 0 ; col < n_col ; col++)
1.1753 + {
1.1754 + last_row = -1 ;
1.1755 +
1.1756 + cp = &A [p [col]] ;
1.1757 + cp_end = &A [p [col+1]] ;
1.1758 +
1.1759 + while (cp < cp_end)
1.1760 + {
1.1761 + row = *cp++ ;
1.1762 +
1.1763 + /* make sure row indices within range */
1.1764 + if (row < 0 || row >= n_row)
1.1765 + {
1.1766 + stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1.1767 + stats [COLAMD_INFO1] = col ;
1.1768 + stats [COLAMD_INFO2] = row ;
1.1769 + stats [COLAMD_INFO3] = n_row ;
1.1770 + DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
1.1771 + return (FALSE) ;
1.1772 + }
1.1773 +
1.1774 + if (row <= last_row || Row [row].shared2.mark == col)
1.1775 + {
1.1776 + /* row index are unsorted or repeated (or both), thus col */
1.1777 + /* is jumbled. This is a notice, not an error condition. */
1.1778 + stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1.1779 + stats [COLAMD_INFO1] = col ;
1.1780 + stats [COLAMD_INFO2] = row ;
1.1781 + (stats [COLAMD_INFO3]) ++ ;
1.1782 + DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
1.1783 + }
1.1784 +
1.1785 + if (Row [row].shared2.mark != col)
1.1786 + {
1.1787 + Row [row].length++ ;
1.1788 + }
1.1789 + else
1.1790 + {
1.1791 + /* this is a repeated entry in the column, */
1.1792 + /* it will be removed */
1.1793 + Col [col].length-- ;
1.1794 + }
1.1795 +
1.1796 + /* mark the row as having been seen in this column */
1.1797 + Row [row].shared2.mark = col ;
1.1798 +
1.1799 + last_row = row ;
1.1800 + }
1.1801 + }
1.1802 +
1.1803 + /* === Compute row pointers ============================================= */
1.1804 +
1.1805 + /* row form of the matrix starts directly after the column */
1.1806 + /* form of matrix in A */
1.1807 + Row [0].start = p [n_col] ;
1.1808 + Row [0].shared1.p = Row [0].start ;
1.1809 + Row [0].shared2.mark = -1 ;
1.1810 + for (row = 1 ; row < n_row ; row++)
1.1811 + {
1.1812 + Row [row].start = Row [row-1].start + Row [row-1].length ;
1.1813 + Row [row].shared1.p = Row [row].start ;
1.1814 + Row [row].shared2.mark = -1 ;
1.1815 + }
1.1816 +
1.1817 + /* === Create row form ================================================== */
1.1818 +
1.1819 + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1.1820 + {
1.1821 + /* if cols jumbled, watch for repeated row indices */
1.1822 + for (col = 0 ; col < n_col ; col++)
1.1823 + {
1.1824 + cp = &A [p [col]] ;
1.1825 + cp_end = &A [p [col+1]] ;
1.1826 + while (cp < cp_end)
1.1827 + {
1.1828 + row = *cp++ ;
1.1829 + if (Row [row].shared2.mark != col)
1.1830 + {
1.1831 + A [(Row [row].shared1.p)++] = col ;
1.1832 + Row [row].shared2.mark = col ;
1.1833 + }
1.1834 + }
1.1835 + }
1.1836 + }
1.1837 + else
1.1838 + {
1.1839 + /* if cols not jumbled, we don't need the mark (this is faster) */
1.1840 + for (col = 0 ; col < n_col ; col++)
1.1841 + {
1.1842 + cp = &A [p [col]] ;
1.1843 + cp_end = &A [p [col+1]] ;
1.1844 + while (cp < cp_end)
1.1845 + {
1.1846 + A [(Row [*cp++].shared1.p)++] = col ;
1.1847 + }
1.1848 + }
1.1849 + }
1.1850 +
1.1851 + /* === Clear the row marks and set row degrees ========================== */
1.1852 +
1.1853 + for (row = 0 ; row < n_row ; row++)
1.1854 + {
1.1855 + Row [row].shared2.mark = 0 ;
1.1856 + Row [row].shared1.degree = Row [row].length ;
1.1857 + }
1.1858 +
1.1859 + /* === See if we need to re-create columns ============================== */
1.1860 +
1.1861 + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1.1862 + {
1.1863 + DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
1.1864 +
1.1865 +#ifndef NDEBUG
1.1866 + /* make sure column lengths are correct */
1.1867 + for (col = 0 ; col < n_col ; col++)
1.1868 + {
1.1869 + p [col] = Col [col].length ;
1.1870 + }
1.1871 + for (row = 0 ; row < n_row ; row++)
1.1872 + {
1.1873 + rp = &A [Row [row].start] ;
1.1874 + rp_end = rp + Row [row].length ;
1.1875 + while (rp < rp_end)
1.1876 + {
1.1877 + p [*rp++]-- ;
1.1878 + }
1.1879 + }
1.1880 + for (col = 0 ; col < n_col ; col++)
1.1881 + {
1.1882 + ASSERT (p [col] == 0) ;
1.1883 + }
1.1884 + /* now p is all zero (different than when debugging is turned off) */
1.1885 +#endif /* NDEBUG */
1.1886 +
1.1887 + /* === Compute col pointers ========================================= */
1.1888 +
1.1889 + /* col form of the matrix starts at A [0]. */
1.1890 + /* Note, we may have a gap between the col form and the row */
1.1891 + /* form if there were duplicate entries, if so, it will be */
1.1892 + /* removed upon the first garbage collection */
1.1893 + Col [0].start = 0 ;
1.1894 + p [0] = Col [0].start ;
1.1895 + for (col = 1 ; col < n_col ; col++)
1.1896 + {
1.1897 + /* note that the lengths here are for pruned columns, i.e. */
1.1898 + /* no duplicate row indices will exist for these columns */
1.1899 + Col [col].start = Col [col-1].start + Col [col-1].length ;
1.1900 + p [col] = Col [col].start ;
1.1901 + }
1.1902 +
1.1903 + /* === Re-create col form =========================================== */
1.1904 +
1.1905 + for (row = 0 ; row < n_row ; row++)
1.1906 + {
1.1907 + rp = &A [Row [row].start] ;
1.1908 + rp_end = rp + Row [row].length ;
1.1909 + while (rp < rp_end)
1.1910 + {
1.1911 + A [(p [*rp++])++] = row ;
1.1912 + }
1.1913 + }
1.1914 + }
1.1915 +
1.1916 + /* === Done. Matrix is not (or no longer) jumbled ====================== */
1.1917 +
1.1918 + return (TRUE) ;
1.1919 +}
1.1920 +
1.1921 +
1.1922 +/* ========================================================================== */
1.1923 +/* === init_scoring ========================================================= */
1.1924 +/* ========================================================================== */
1.1925 +
1.1926 +/*
1.1927 + Kills dense or empty columns and rows, calculates an initial score for
1.1928 + each column, and places all columns in the degree lists. Not user-callable.
1.1929 +*/
1.1930 +
1.1931 +PRIVATE void init_scoring
1.1932 +(
1.1933 + /* === Parameters ======================================================= */
1.1934 +
1.1935 + Int n_row, /* number of rows of A */
1.1936 + Int n_col, /* number of columns of A */
1.1937 + Colamd_Row Row [], /* of size n_row+1 */
1.1938 + Colamd_Col Col [], /* of size n_col+1 */
1.1939 + Int A [], /* column form and row form of A */
1.1940 + Int head [], /* of size n_col+1 */
1.1941 + double knobs [COLAMD_KNOBS],/* parameters */
1.1942 + Int *p_n_row2, /* number of non-dense, non-empty rows */
1.1943 + Int *p_n_col2, /* number of non-dense, non-empty columns */
1.1944 + Int *p_max_deg /* maximum row degree */
1.1945 +)
1.1946 +{
1.1947 + /* === Local variables ================================================== */
1.1948 +
1.1949 + Int c ; /* a column index */
1.1950 + Int r, row ; /* a row index */
1.1951 + Int *cp ; /* a column pointer */
1.1952 + Int deg ; /* degree of a row or column */
1.1953 + Int *cp_end ; /* a pointer to the end of a column */
1.1954 + Int *new_cp ; /* new column pointer */
1.1955 + Int col_length ; /* length of pruned column */
1.1956 + Int score ; /* current column score */
1.1957 + Int n_col2 ; /* number of non-dense, non-empty columns */
1.1958 + Int n_row2 ; /* number of non-dense, non-empty rows */
1.1959 + Int dense_row_count ; /* remove rows with more entries than this */
1.1960 + Int dense_col_count ; /* remove cols with more entries than this */
1.1961 + Int min_score ; /* smallest column score */
1.1962 + Int max_deg ; /* maximum row degree */
1.1963 + Int next_col ; /* Used to add to degree list.*/
1.1964 +
1.1965 +#ifndef NDEBUG
1.1966 + Int debug_count ; /* debug only. */
1.1967 +#endif /* NDEBUG */
1.1968 +
1.1969 + /* === Extract knobs ==================================================== */
1.1970 +
1.1971 + /* Note: if knobs contains a NaN, this is undefined: */
1.1972 + if (knobs [COLAMD_DENSE_ROW] < 0)
1.1973 + {
1.1974 + /* only remove completely dense rows */
1.1975 + dense_row_count = n_col-1 ;
1.1976 + }
1.1977 + else
1.1978 + {
1.1979 + dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
1.1980 + }
1.1981 + if (knobs [COLAMD_DENSE_COL] < 0)
1.1982 + {
1.1983 + /* only remove completely dense columns */
1.1984 + dense_col_count = n_row-1 ;
1.1985 + }
1.1986 + else
1.1987 + {
1.1988 + dense_col_count =
1.1989 + DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
1.1990 + }
1.1991 +
1.1992 + DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
1.1993 + max_deg = 0 ;
1.1994 + n_col2 = n_col ;
1.1995 + n_row2 = n_row ;
1.1996 +
1.1997 + /* === Kill empty columns =============================================== */
1.1998 +
1.1999 + /* Put the empty columns at the end in their natural order, so that LU */
1.2000 + /* factorization can proceed as far as possible. */
1.2001 + for (c = n_col-1 ; c >= 0 ; c--)
1.2002 + {
1.2003 + deg = Col [c].length ;
1.2004 + if (deg == 0)
1.2005 + {
1.2006 + /* this is a empty column, kill and order it last */
1.2007 + Col [c].shared2.order = --n_col2 ;
1.2008 + KILL_PRINCIPAL_COL (c) ;
1.2009 + }
1.2010 + }
1.2011 + DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
1.2012 +
1.2013 + /* === Kill dense columns =============================================== */
1.2014 +
1.2015 + /* Put the dense columns at the end, in their natural order */
1.2016 + for (c = n_col-1 ; c >= 0 ; c--)
1.2017 + {
1.2018 + /* skip any dead columns */
1.2019 + if (COL_IS_DEAD (c))
1.2020 + {
1.2021 + continue ;
1.2022 + }
1.2023 + deg = Col [c].length ;
1.2024 + if (deg > dense_col_count)
1.2025 + {
1.2026 + /* this is a dense column, kill and order it last */
1.2027 + Col [c].shared2.order = --n_col2 ;
1.2028 + /* decrement the row degrees */
1.2029 + cp = &A [Col [c].start] ;
1.2030 + cp_end = cp + Col [c].length ;
1.2031 + while (cp < cp_end)
1.2032 + {
1.2033 + Row [*cp++].shared1.degree-- ;
1.2034 + }
1.2035 + KILL_PRINCIPAL_COL (c) ;
1.2036 + }
1.2037 + }
1.2038 + DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
1.2039 +
1.2040 + /* === Kill dense and empty rows ======================================== */
1.2041 +
1.2042 + for (r = 0 ; r < n_row ; r++)
1.2043 + {
1.2044 + deg = Row [r].shared1.degree ;
1.2045 + ASSERT (deg >= 0 && deg <= n_col) ;
1.2046 + if (deg > dense_row_count || deg == 0)
1.2047 + {
1.2048 + /* kill a dense or empty row */
1.2049 + KILL_ROW (r) ;
1.2050 + --n_row2 ;
1.2051 + }
1.2052 + else
1.2053 + {
1.2054 + /* keep track of max degree of remaining rows */
1.2055 + max_deg = MAX (max_deg, deg) ;
1.2056 + }
1.2057 + }
1.2058 + DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
1.2059 +
1.2060 + /* === Compute initial column scores ==================================== */
1.2061 +
1.2062 + /* At this point the row degrees are accurate. They reflect the number */
1.2063 + /* of "live" (non-dense) columns in each row. No empty rows exist. */
1.2064 + /* Some "live" columns may contain only dead rows, however. These are */
1.2065 + /* pruned in the code below. */
1.2066 +
1.2067 + /* now find the initial matlab score for each column */
1.2068 + for (c = n_col-1 ; c >= 0 ; c--)
1.2069 + {
1.2070 + /* skip dead column */
1.2071 + if (COL_IS_DEAD (c))
1.2072 + {
1.2073 + continue ;
1.2074 + }
1.2075 + score = 0 ;
1.2076 + cp = &A [Col [c].start] ;
1.2077 + new_cp = cp ;
1.2078 + cp_end = cp + Col [c].length ;
1.2079 + while (cp < cp_end)
1.2080 + {
1.2081 + /* get a row */
1.2082 + row = *cp++ ;
1.2083 + /* skip if dead */
1.2084 + if (ROW_IS_DEAD (row))
1.2085 + {
1.2086 + continue ;
1.2087 + }
1.2088 + /* compact the column */
1.2089 + *new_cp++ = row ;
1.2090 + /* add row's external degree */
1.2091 + score += Row [row].shared1.degree - 1 ;
1.2092 + /* guard against integer overflow */
1.2093 + score = MIN (score, n_col) ;
1.2094 + }
1.2095 + /* determine pruned column length */
1.2096 + col_length = (Int) (new_cp - &A [Col [c].start]) ;
1.2097 + if (col_length == 0)
1.2098 + {
1.2099 + /* a newly-made null column (all rows in this col are "dense" */
1.2100 + /* and have already been killed) */
1.2101 + DEBUG2 (("Newly null killed: %d\n", c)) ;
1.2102 + Col [c].shared2.order = --n_col2 ;
1.2103 + KILL_PRINCIPAL_COL (c) ;
1.2104 + }
1.2105 + else
1.2106 + {
1.2107 + /* set column length and set score */
1.2108 + ASSERT (score >= 0) ;
1.2109 + ASSERT (score <= n_col) ;
1.2110 + Col [c].length = col_length ;
1.2111 + Col [c].shared2.score = score ;
1.2112 + }
1.2113 + }
1.2114 + DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
1.2115 + n_col-n_col2)) ;
1.2116 +
1.2117 + /* At this point, all empty rows and columns are dead. All live columns */
1.2118 + /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
1.2119 + /* yet). Rows may contain dead columns, but all live rows contain at */
1.2120 + /* least one live column. */
1.2121 +
1.2122 +#ifndef NDEBUG
1.2123 + debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
1.2124 +#endif /* NDEBUG */
1.2125 +
1.2126 + /* === Initialize degree lists ========================================== */
1.2127 +
1.2128 +#ifndef NDEBUG
1.2129 + debug_count = 0 ;
1.2130 +#endif /* NDEBUG */
1.2131 +
1.2132 + /* clear the hash buckets */
1.2133 + for (c = 0 ; c <= n_col ; c++)
1.2134 + {
1.2135 + head [c] = EMPTY ;
1.2136 + }
1.2137 + min_score = n_col ;
1.2138 + /* place in reverse order, so low column indices are at the front */
1.2139 + /* of the lists. This is to encourage natural tie-breaking */
1.2140 + for (c = n_col-1 ; c >= 0 ; c--)
1.2141 + {
1.2142 + /* only add principal columns to degree lists */
1.2143 + if (COL_IS_ALIVE (c))
1.2144 + {
1.2145 + DEBUG4 (("place %d score %d minscore %d ncol %d\n",
1.2146 + c, Col [c].shared2.score, min_score, n_col)) ;
1.2147 +
1.2148 + /* === Add columns score to DList =============================== */
1.2149 +
1.2150 + score = Col [c].shared2.score ;
1.2151 +
1.2152 + ASSERT (min_score >= 0) ;
1.2153 + ASSERT (min_score <= n_col) ;
1.2154 + ASSERT (score >= 0) ;
1.2155 + ASSERT (score <= n_col) ;
1.2156 + ASSERT (head [score] >= EMPTY) ;
1.2157 +
1.2158 + /* now add this column to dList at proper score location */
1.2159 + next_col = head [score] ;
1.2160 + Col [c].shared3.prev = EMPTY ;
1.2161 + Col [c].shared4.degree_next = next_col ;
1.2162 +
1.2163 + /* if there already was a column with the same score, set its */
1.2164 + /* previous pointer to this new column */
1.2165 + if (next_col != EMPTY)
1.2166 + {
1.2167 + Col [next_col].shared3.prev = c ;
1.2168 + }
1.2169 + head [score] = c ;
1.2170 +
1.2171 + /* see if this score is less than current min */
1.2172 + min_score = MIN (min_score, score) ;
1.2173 +
1.2174 +#ifndef NDEBUG
1.2175 + debug_count++ ;
1.2176 +#endif /* NDEBUG */
1.2177 +
1.2178 + }
1.2179 + }
1.2180 +
1.2181 +#ifndef NDEBUG
1.2182 + DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
1.2183 + debug_count, n_col, n_col-debug_count)) ;
1.2184 + ASSERT (debug_count == n_col2) ;
1.2185 + debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
1.2186 +#endif /* NDEBUG */
1.2187 +
1.2188 + /* === Return number of remaining columns, and max row degree =========== */
1.2189 +
1.2190 + *p_n_col2 = n_col2 ;
1.2191 + *p_n_row2 = n_row2 ;
1.2192 + *p_max_deg = max_deg ;
1.2193 +}
1.2194 +
1.2195 +
1.2196 +/* ========================================================================== */
1.2197 +/* === find_ordering ======================================================== */
1.2198 +/* ========================================================================== */
1.2199 +
1.2200 +/*
1.2201 + Order the principal columns of the supercolumn form of the matrix
1.2202 + (no supercolumns on input). Uses a minimum approximate column minimum
1.2203 + degree ordering method. Not user-callable.
1.2204 +*/
1.2205 +
1.2206 +PRIVATE Int find_ordering /* return the number of garbage collections */
1.2207 +(
1.2208 + /* === Parameters ======================================================= */
1.2209 +
1.2210 + Int n_row, /* number of rows of A */
1.2211 + Int n_col, /* number of columns of A */
1.2212 + Int Alen, /* size of A, 2*nnz + n_col or larger */
1.2213 + Colamd_Row Row [], /* of size n_row+1 */
1.2214 + Colamd_Col Col [], /* of size n_col+1 */
1.2215 + Int A [], /* column form and row form of A */
1.2216 + Int head [], /* of size n_col+1 */
1.2217 + Int n_col2, /* Remaining columns to order */
1.2218 + Int max_deg, /* Maximum row degree */
1.2219 + Int pfree, /* index of first free slot (2*nnz on entry) */
1.2220 + Int aggressive
1.2221 +)
1.2222 +{
1.2223 + /* === Local variables ================================================== */
1.2224 +
1.2225 + Int k ; /* current pivot ordering step */
1.2226 + Int pivot_col ; /* current pivot column */
1.2227 + Int *cp ; /* a column pointer */
1.2228 + Int *rp ; /* a row pointer */
1.2229 + Int pivot_row ; /* current pivot row */
1.2230 + Int *new_cp ; /* modified column pointer */
1.2231 + Int *new_rp ; /* modified row pointer */
1.2232 + Int pivot_row_start ; /* pointer to start of pivot row */
1.2233 + Int pivot_row_degree ; /* number of columns in pivot row */
1.2234 + Int pivot_row_length ; /* number of supercolumns in pivot row */
1.2235 + Int pivot_col_score ; /* score of pivot column */
1.2236 + Int needed_memory ; /* free space needed for pivot row */
1.2237 + Int *cp_end ; /* pointer to the end of a column */
1.2238 + Int *rp_end ; /* pointer to the end of a row */
1.2239 + Int row ; /* a row index */
1.2240 + Int col ; /* a column index */
1.2241 + Int max_score ; /* maximum possible score */
1.2242 + Int cur_score ; /* score of current column */
1.2243 + unsigned Int hash ; /* hash value for supernode detection */
1.2244 + Int head_column ; /* head of hash bucket */
1.2245 + Int first_col ; /* first column in hash bucket */
1.2246 + Int tag_mark ; /* marker value for mark array */
1.2247 + Int row_mark ; /* Row [row].shared2.mark */
1.2248 + Int set_difference ; /* set difference size of row with pivot row */
1.2249 + Int min_score ; /* smallest column score */
1.2250 + Int col_thickness ; /* "thickness" (no. of columns in a supercol) */
1.2251 + Int max_mark ; /* maximum value of tag_mark */
1.2252 + Int pivot_col_thickness ; /* number of columns represented by pivot col */
1.2253 + Int prev_col ; /* Used by Dlist operations. */
1.2254 + Int next_col ; /* Used by Dlist operations. */
1.2255 + Int ngarbage ; /* number of garbage collections performed */
1.2256 +
1.2257 +#ifndef NDEBUG
1.2258 + Int debug_d ; /* debug loop counter */
1.2259 + Int debug_step = 0 ; /* debug loop counter */
1.2260 +#endif /* NDEBUG */
1.2261 +
1.2262 + /* === Initialization and clear mark ==================================== */
1.2263 +
1.2264 + max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */
1.2265 + tag_mark = clear_mark (0, max_mark, n_row, Row) ;
1.2266 + min_score = 0 ;
1.2267 + ngarbage = 0 ;
1.2268 + DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
1.2269 +
1.2270 + /* === Order the columns ================================================ */
1.2271 +
1.2272 + for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
1.2273 + {
1.2274 +
1.2275 +#ifndef NDEBUG
1.2276 + if (debug_step % 100 == 0)
1.2277 + {
1.2278 + DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ;
1.2279 + }
1.2280 + else
1.2281 + {
1.2282 + DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
1.2283 + }
1.2284 + debug_step++ ;
1.2285 + debug_deg_lists (n_row, n_col, Row, Col, head,
1.2286 + min_score, n_col2-k, max_deg) ;
1.2287 + debug_matrix (n_row, n_col, Row, Col, A) ;
1.2288 +#endif /* NDEBUG */
1.2289 +
1.2290 + /* === Select pivot column, and order it ============================ */
1.2291 +
1.2292 + /* make sure degree list isn't empty */
1.2293 + ASSERT (min_score >= 0) ;
1.2294 + ASSERT (min_score <= n_col) ;
1.2295 + ASSERT (head [min_score] >= EMPTY) ;
1.2296 +
1.2297 +#ifndef NDEBUG
1.2298 + for (debug_d = 0 ; debug_d < min_score ; debug_d++)
1.2299 + {
1.2300 + ASSERT (head [debug_d] == EMPTY) ;
1.2301 + }
1.2302 +#endif /* NDEBUG */
1.2303 +
1.2304 + /* get pivot column from head of minimum degree list */
1.2305 + while (head [min_score] == EMPTY && min_score < n_col)
1.2306 + {
1.2307 + min_score++ ;
1.2308 + }
1.2309 + pivot_col = head [min_score] ;
1.2310 + ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
1.2311 + next_col = Col [pivot_col].shared4.degree_next ;
1.2312 + head [min_score] = next_col ;
1.2313 + if (next_col != EMPTY)
1.2314 + {
1.2315 + Col [next_col].shared3.prev = EMPTY ;
1.2316 + }
1.2317 +
1.2318 + ASSERT (COL_IS_ALIVE (pivot_col)) ;
1.2319 +
1.2320 + /* remember score for defrag check */
1.2321 + pivot_col_score = Col [pivot_col].shared2.score ;
1.2322 +
1.2323 + /* the pivot column is the kth column in the pivot order */
1.2324 + Col [pivot_col].shared2.order = k ;
1.2325 +
1.2326 + /* increment order count by column thickness */
1.2327 + pivot_col_thickness = Col [pivot_col].shared1.thickness ;
1.2328 + k += pivot_col_thickness ;
1.2329 + ASSERT (pivot_col_thickness > 0) ;
1.2330 + DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
1.2331 +
1.2332 + /* === Garbage_collection, if necessary ============================= */
1.2333 +
1.2334 + needed_memory = MIN (pivot_col_score, n_col - k) ;
1.2335 + if (pfree + needed_memory >= Alen)
1.2336 + {
1.2337 + pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
1.2338 + ngarbage++ ;
1.2339 + /* after garbage collection we will have enough */
1.2340 + ASSERT (pfree + needed_memory < Alen) ;
1.2341 + /* garbage collection has wiped out the Row[].shared2.mark array */
1.2342 + tag_mark = clear_mark (0, max_mark, n_row, Row) ;
1.2343 +
1.2344 +#ifndef NDEBUG
1.2345 + debug_matrix (n_row, n_col, Row, Col, A) ;
1.2346 +#endif /* NDEBUG */
1.2347 + }
1.2348 +
1.2349 + /* === Compute pivot row pattern ==================================== */
1.2350 +
1.2351 + /* get starting location for this new merged row */
1.2352 + pivot_row_start = pfree ;
1.2353 +
1.2354 + /* initialize new row counts to zero */
1.2355 + pivot_row_degree = 0 ;
1.2356 +
1.2357 + /* tag pivot column as having been visited so it isn't included */
1.2358 + /* in merged pivot row */
1.2359 + Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
1.2360 +
1.2361 + /* pivot row is the union of all rows in the pivot column pattern */
1.2362 + cp = &A [Col [pivot_col].start] ;
1.2363 + cp_end = cp + Col [pivot_col].length ;
1.2364 + while (cp < cp_end)
1.2365 + {
1.2366 + /* get a row */
1.2367 + row = *cp++ ;
1.2368 + DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
1.2369 + /* skip if row is dead */
1.2370 + if (ROW_IS_ALIVE (row))
1.2371 + {
1.2372 + rp = &A [Row [row].start] ;
1.2373 + rp_end = rp + Row [row].length ;
1.2374 + while (rp < rp_end)
1.2375 + {
1.2376 + /* get a column */
1.2377 + col = *rp++ ;
1.2378 + /* add the column, if alive and untagged */
1.2379 + col_thickness = Col [col].shared1.thickness ;
1.2380 + if (col_thickness > 0 && COL_IS_ALIVE (col))
1.2381 + {
1.2382 + /* tag column in pivot row */
1.2383 + Col [col].shared1.thickness = -col_thickness ;
1.2384 + ASSERT (pfree < Alen) ;
1.2385 + /* place column in pivot row */
1.2386 + A [pfree++] = col ;
1.2387 + pivot_row_degree += col_thickness ;
1.2388 + }
1.2389 + }
1.2390 + }
1.2391 + }
1.2392 +
1.2393 + /* clear tag on pivot column */
1.2394 + Col [pivot_col].shared1.thickness = pivot_col_thickness ;
1.2395 + max_deg = MAX (max_deg, pivot_row_degree) ;
1.2396 +
1.2397 +#ifndef NDEBUG
1.2398 + DEBUG3 (("check2\n")) ;
1.2399 + debug_mark (n_row, Row, tag_mark, max_mark) ;
1.2400 +#endif /* NDEBUG */
1.2401 +
1.2402 + /* === Kill all rows used to construct pivot row ==================== */
1.2403 +
1.2404 + /* also kill pivot row, temporarily */
1.2405 + cp = &A [Col [pivot_col].start] ;
1.2406 + cp_end = cp + Col [pivot_col].length ;
1.2407 + while (cp < cp_end)
1.2408 + {
1.2409 + /* may be killing an already dead row */
1.2410 + row = *cp++ ;
1.2411 + DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
1.2412 + KILL_ROW (row) ;
1.2413 + }
1.2414 +
1.2415 + /* === Select a row index to use as the new pivot row =============== */
1.2416 +
1.2417 + pivot_row_length = pfree - pivot_row_start ;
1.2418 + if (pivot_row_length > 0)
1.2419 + {
1.2420 + /* pick the "pivot" row arbitrarily (first row in col) */
1.2421 + pivot_row = A [Col [pivot_col].start] ;
1.2422 + DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
1.2423 + }
1.2424 + else
1.2425 + {
1.2426 + /* there is no pivot row, since it is of zero length */
1.2427 + pivot_row = EMPTY ;
1.2428 + ASSERT (pivot_row_length == 0) ;
1.2429 + }
1.2430 + ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
1.2431 +
1.2432 + /* === Approximate degree computation =============================== */
1.2433 +
1.2434 + /* Here begins the computation of the approximate degree. The column */
1.2435 + /* score is the sum of the pivot row "length", plus the size of the */
1.2436 + /* set differences of each row in the column minus the pattern of the */
1.2437 + /* pivot row itself. The column ("thickness") itself is also */
1.2438 + /* excluded from the column score (we thus use an approximate */
1.2439 + /* external degree). */
1.2440 +
1.2441 + /* The time taken by the following code (compute set differences, and */
1.2442 + /* add them up) is proportional to the size of the data structure */
1.2443 + /* being scanned - that is, the sum of the sizes of each column in */
1.2444 + /* the pivot row. Thus, the amortized time to compute a column score */
1.2445 + /* is proportional to the size of that column (where size, in this */
1.2446 + /* context, is the column "length", or the number of row indices */
1.2447 + /* in that column). The number of row indices in a column is */
1.2448 + /* monotonically non-decreasing, from the length of the original */
1.2449 + /* column on input to colamd. */
1.2450 +
1.2451 + /* === Compute set differences ====================================== */
1.2452 +
1.2453 + DEBUG3 (("** Computing set differences phase. **\n")) ;
1.2454 +
1.2455 + /* pivot row is currently dead - it will be revived later. */
1.2456 +
1.2457 + DEBUG3 (("Pivot row: ")) ;
1.2458 + /* for each column in pivot row */
1.2459 + rp = &A [pivot_row_start] ;
1.2460 + rp_end = rp + pivot_row_length ;
1.2461 + while (rp < rp_end)
1.2462 + {
1.2463 + col = *rp++ ;
1.2464 + ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1.2465 + DEBUG3 (("Col: %d\n", col)) ;
1.2466 +
1.2467 + /* clear tags used to construct pivot row pattern */
1.2468 + col_thickness = -Col [col].shared1.thickness ;
1.2469 + ASSERT (col_thickness > 0) ;
1.2470 + Col [col].shared1.thickness = col_thickness ;
1.2471 +
1.2472 + /* === Remove column from degree list =========================== */
1.2473 +
1.2474 + cur_score = Col [col].shared2.score ;
1.2475 + prev_col = Col [col].shared3.prev ;
1.2476 + next_col = Col [col].shared4.degree_next ;
1.2477 + ASSERT (cur_score >= 0) ;
1.2478 + ASSERT (cur_score <= n_col) ;
1.2479 + ASSERT (cur_score >= EMPTY) ;
1.2480 + if (prev_col == EMPTY)
1.2481 + {
1.2482 + head [cur_score] = next_col ;
1.2483 + }
1.2484 + else
1.2485 + {
1.2486 + Col [prev_col].shared4.degree_next = next_col ;
1.2487 + }
1.2488 + if (next_col != EMPTY)
1.2489 + {
1.2490 + Col [next_col].shared3.prev = prev_col ;
1.2491 + }
1.2492 +
1.2493 + /* === Scan the column ========================================== */
1.2494 +
1.2495 + cp = &A [Col [col].start] ;
1.2496 + cp_end = cp + Col [col].length ;
1.2497 + while (cp < cp_end)
1.2498 + {
1.2499 + /* get a row */
1.2500 + row = *cp++ ;
1.2501 + row_mark = Row [row].shared2.mark ;
1.2502 + /* skip if dead */
1.2503 + if (ROW_IS_MARKED_DEAD (row_mark))
1.2504 + {
1.2505 + continue ;
1.2506 + }
1.2507 + ASSERT (row != pivot_row) ;
1.2508 + set_difference = row_mark - tag_mark ;
1.2509 + /* check if the row has been seen yet */
1.2510 + if (set_difference < 0)
1.2511 + {
1.2512 + ASSERT (Row [row].shared1.degree <= max_deg) ;
1.2513 + set_difference = Row [row].shared1.degree ;
1.2514 + }
1.2515 + /* subtract column thickness from this row's set difference */
1.2516 + set_difference -= col_thickness ;
1.2517 + ASSERT (set_difference >= 0) ;
1.2518 + /* absorb this row if the set difference becomes zero */
1.2519 + if (set_difference == 0 && aggressive)
1.2520 + {
1.2521 + DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
1.2522 + KILL_ROW (row) ;
1.2523 + }
1.2524 + else
1.2525 + {
1.2526 + /* save the new mark */
1.2527 + Row [row].shared2.mark = set_difference + tag_mark ;
1.2528 + }
1.2529 + }
1.2530 + }
1.2531 +
1.2532 +#ifndef NDEBUG
1.2533 + debug_deg_lists (n_row, n_col, Row, Col, head,
1.2534 + min_score, n_col2-k-pivot_row_degree, max_deg) ;
1.2535 +#endif /* NDEBUG */
1.2536 +
1.2537 + /* === Add up set differences for each column ======================= */
1.2538 +
1.2539 + DEBUG3 (("** Adding set differences phase. **\n")) ;
1.2540 +
1.2541 + /* for each column in pivot row */
1.2542 + rp = &A [pivot_row_start] ;
1.2543 + rp_end = rp + pivot_row_length ;
1.2544 + while (rp < rp_end)
1.2545 + {
1.2546 + /* get a column */
1.2547 + col = *rp++ ;
1.2548 + ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1.2549 + hash = 0 ;
1.2550 + cur_score = 0 ;
1.2551 + cp = &A [Col [col].start] ;
1.2552 + /* compact the column */
1.2553 + new_cp = cp ;
1.2554 + cp_end = cp + Col [col].length ;
1.2555 +
1.2556 + DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
1.2557 +
1.2558 + while (cp < cp_end)
1.2559 + {
1.2560 + /* get a row */
1.2561 + row = *cp++ ;
1.2562 + ASSERT(row >= 0 && row < n_row) ;
1.2563 + row_mark = Row [row].shared2.mark ;
1.2564 + /* skip if dead */
1.2565 + if (ROW_IS_MARKED_DEAD (row_mark))
1.2566 + {
1.2567 + DEBUG4 ((" Row %d, dead\n", row)) ;
1.2568 + continue ;
1.2569 + }
1.2570 + DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
1.2571 + ASSERT (row_mark >= tag_mark) ;
1.2572 + /* compact the column */
1.2573 + *new_cp++ = row ;
1.2574 + /* compute hash function */
1.2575 + hash += row ;
1.2576 + /* add set difference */
1.2577 + cur_score += row_mark - tag_mark ;
1.2578 + /* integer overflow... */
1.2579 + cur_score = MIN (cur_score, n_col) ;
1.2580 + }
1.2581 +
1.2582 + /* recompute the column's length */
1.2583 + Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
1.2584 +
1.2585 + /* === Further mass elimination ================================= */
1.2586 +
1.2587 + if (Col [col].length == 0)
1.2588 + {
1.2589 + DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
1.2590 + /* nothing left but the pivot row in this column */
1.2591 + KILL_PRINCIPAL_COL (col) ;
1.2592 + pivot_row_degree -= Col [col].shared1.thickness ;
1.2593 + ASSERT (pivot_row_degree >= 0) ;
1.2594 + /* order it */
1.2595 + Col [col].shared2.order = k ;
1.2596 + /* increment order count by column thickness */
1.2597 + k += Col [col].shared1.thickness ;
1.2598 + }
1.2599 + else
1.2600 + {
1.2601 + /* === Prepare for supercolumn detection ==================== */
1.2602 +
1.2603 + DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
1.2604 +
1.2605 + /* save score so far */
1.2606 + Col [col].shared2.score = cur_score ;
1.2607 +
1.2608 + /* add column to hash table, for supercolumn detection */
1.2609 + hash %= n_col + 1 ;
1.2610 +
1.2611 + DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
1.2612 + ASSERT (((Int) hash) <= n_col) ;
1.2613 +
1.2614 + head_column = head [hash] ;
1.2615 + if (head_column > EMPTY)
1.2616 + {
1.2617 + /* degree list "hash" is non-empty, use prev (shared3) of */
1.2618 + /* first column in degree list as head of hash bucket */
1.2619 + first_col = Col [head_column].shared3.headhash ;
1.2620 + Col [head_column].shared3.headhash = col ;
1.2621 + }
1.2622 + else
1.2623 + {
1.2624 + /* degree list "hash" is empty, use head as hash bucket */
1.2625 + first_col = - (head_column + 2) ;
1.2626 + head [hash] = - (col + 2) ;
1.2627 + }
1.2628 + Col [col].shared4.hash_next = first_col ;
1.2629 +
1.2630 + /* save hash function in Col [col].shared3.hash */
1.2631 + Col [col].shared3.hash = (Int) hash ;
1.2632 + ASSERT (COL_IS_ALIVE (col)) ;
1.2633 + }
1.2634 + }
1.2635 +
1.2636 + /* The approximate external column degree is now computed. */
1.2637 +
1.2638 + /* === Supercolumn detection ======================================== */
1.2639 +
1.2640 + DEBUG3 (("** Supercolumn detection phase. **\n")) ;
1.2641 +
1.2642 + detect_super_cols (
1.2643 +
1.2644 +#ifndef NDEBUG
1.2645 + n_col, Row,
1.2646 +#endif /* NDEBUG */
1.2647 +
1.2648 + Col, A, head, pivot_row_start, pivot_row_length) ;
1.2649 +
1.2650 + /* === Kill the pivotal column ====================================== */
1.2651 +
1.2652 + KILL_PRINCIPAL_COL (pivot_col) ;
1.2653 +
1.2654 + /* === Clear mark =================================================== */
1.2655 +
1.2656 + tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
1.2657 +
1.2658 +#ifndef NDEBUG
1.2659 + DEBUG3 (("check3\n")) ;
1.2660 + debug_mark (n_row, Row, tag_mark, max_mark) ;
1.2661 +#endif /* NDEBUG */
1.2662 +
1.2663 + /* === Finalize the new pivot row, and column scores ================ */
1.2664 +
1.2665 + DEBUG3 (("** Finalize scores phase. **\n")) ;
1.2666 +
1.2667 + /* for each column in pivot row */
1.2668 + rp = &A [pivot_row_start] ;
1.2669 + /* compact the pivot row */
1.2670 + new_rp = rp ;
1.2671 + rp_end = rp + pivot_row_length ;
1.2672 + while (rp < rp_end)
1.2673 + {
1.2674 + col = *rp++ ;
1.2675 + /* skip dead columns */
1.2676 + if (COL_IS_DEAD (col))
1.2677 + {
1.2678 + continue ;
1.2679 + }
1.2680 + *new_rp++ = col ;
1.2681 + /* add new pivot row to column */
1.2682 + A [Col [col].start + (Col [col].length++)] = pivot_row ;
1.2683 +
1.2684 + /* retrieve score so far and add on pivot row's degree. */
1.2685 + /* (we wait until here for this in case the pivot */
1.2686 + /* row's degree was reduced due to mass elimination). */
1.2687 + cur_score = Col [col].shared2.score + pivot_row_degree ;
1.2688 +
1.2689 + /* calculate the max possible score as the number of */
1.2690 + /* external columns minus the 'k' value minus the */
1.2691 + /* columns thickness */
1.2692 + max_score = n_col - k - Col [col].shared1.thickness ;
1.2693 +
1.2694 + /* make the score the external degree of the union-of-rows */
1.2695 + cur_score -= Col [col].shared1.thickness ;
1.2696 +
1.2697 + /* make sure score is less or equal than the max score */
1.2698 + cur_score = MIN (cur_score, max_score) ;
1.2699 + ASSERT (cur_score >= 0) ;
1.2700 +
1.2701 + /* store updated score */
1.2702 + Col [col].shared2.score = cur_score ;
1.2703 +
1.2704 + /* === Place column back in degree list ========================= */
1.2705 +
1.2706 + ASSERT (min_score >= 0) ;
1.2707 + ASSERT (min_score <= n_col) ;
1.2708 + ASSERT (cur_score >= 0) ;
1.2709 + ASSERT (cur_score <= n_col) ;
1.2710 + ASSERT (head [cur_score] >= EMPTY) ;
1.2711 + next_col = head [cur_score] ;
1.2712 + Col [col].shared4.degree_next = next_col ;
1.2713 + Col [col].shared3.prev = EMPTY ;
1.2714 + if (next_col != EMPTY)
1.2715 + {
1.2716 + Col [next_col].shared3.prev = col ;
1.2717 + }
1.2718 + head [cur_score] = col ;
1.2719 +
1.2720 + /* see if this score is less than current min */
1.2721 + min_score = MIN (min_score, cur_score) ;
1.2722 +
1.2723 + }
1.2724 +
1.2725 +#ifndef NDEBUG
1.2726 + debug_deg_lists (n_row, n_col, Row, Col, head,
1.2727 + min_score, n_col2-k, max_deg) ;
1.2728 +#endif /* NDEBUG */
1.2729 +
1.2730 + /* === Resurrect the new pivot row ================================== */
1.2731 +
1.2732 + if (pivot_row_degree > 0)
1.2733 + {
1.2734 + /* update pivot row length to reflect any cols that were killed */
1.2735 + /* during super-col detection and mass elimination */
1.2736 + Row [pivot_row].start = pivot_row_start ;
1.2737 + Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
1.2738 + ASSERT (Row [pivot_row].length > 0) ;
1.2739 + Row [pivot_row].shared1.degree = pivot_row_degree ;
1.2740 + Row [pivot_row].shared2.mark = 0 ;
1.2741 + /* pivot row is no longer dead */
1.2742 +
1.2743 + DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
1.2744 + pivot_row, pivot_row_degree)) ;
1.2745 + }
1.2746 + }
1.2747 +
1.2748 + /* === All principal columns have now been ordered ====================== */
1.2749 +
1.2750 + return (ngarbage) ;
1.2751 +}
1.2752 +
1.2753 +
1.2754 +/* ========================================================================== */
1.2755 +/* === order_children ======================================================= */
1.2756 +/* ========================================================================== */
1.2757 +
1.2758 +/*
1.2759 + The find_ordering routine has ordered all of the principal columns (the
1.2760 + representatives of the supercolumns). The non-principal columns have not
1.2761 + yet been ordered. This routine orders those columns by walking up the
1.2762 + parent tree (a column is a child of the column which absorbed it). The
1.2763 + final permutation vector is then placed in p [0 ... n_col-1], with p [0]
1.2764 + being the first column, and p [n_col-1] being the last. It doesn't look
1.2765 + like it at first glance, but be assured that this routine takes time linear
1.2766 + in the number of columns. Although not immediately obvious, the time
1.2767 + taken by this routine is O (n_col), that is, linear in the number of
1.2768 + columns. Not user-callable.
1.2769 +*/
1.2770 +
1.2771 +PRIVATE void order_children
1.2772 +(
1.2773 + /* === Parameters ======================================================= */
1.2774 +
1.2775 + Int n_col, /* number of columns of A */
1.2776 + Colamd_Col Col [], /* of size n_col+1 */
1.2777 + Int p [] /* p [0 ... n_col-1] is the column permutation*/
1.2778 +)
1.2779 +{
1.2780 + /* === Local variables ================================================== */
1.2781 +
1.2782 + Int i ; /* loop counter for all columns */
1.2783 + Int c ; /* column index */
1.2784 + Int parent ; /* index of column's parent */
1.2785 + Int order ; /* column's order */
1.2786 +
1.2787 + /* === Order each non-principal column ================================== */
1.2788 +
1.2789 + for (i = 0 ; i < n_col ; i++)
1.2790 + {
1.2791 + /* find an un-ordered non-principal column */
1.2792 + ASSERT (COL_IS_DEAD (i)) ;
1.2793 + if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
1.2794 + {
1.2795 + parent = i ;
1.2796 + /* once found, find its principal parent */
1.2797 + do
1.2798 + {
1.2799 + parent = Col [parent].shared1.parent ;
1.2800 + } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
1.2801 +
1.2802 + /* now, order all un-ordered non-principal columns along path */
1.2803 + /* to this parent. collapse tree at the same time */
1.2804 + c = i ;
1.2805 + /* get order of parent */
1.2806 + order = Col [parent].shared2.order ;
1.2807 +
1.2808 + do
1.2809 + {
1.2810 + ASSERT (Col [c].shared2.order == EMPTY) ;
1.2811 +
1.2812 + /* order this column */
1.2813 + Col [c].shared2.order = order++ ;
1.2814 + /* collaps tree */
1.2815 + Col [c].shared1.parent = parent ;
1.2816 +
1.2817 + /* get immediate parent of this column */
1.2818 + c = Col [c].shared1.parent ;
1.2819 +
1.2820 + /* continue until we hit an ordered column. There are */
1.2821 + /* guarranteed not to be anymore unordered columns */
1.2822 + /* above an ordered column */
1.2823 + } while (Col [c].shared2.order == EMPTY) ;
1.2824 +
1.2825 + /* re-order the super_col parent to largest order for this group */
1.2826 + Col [parent].shared2.order = order ;
1.2827 + }
1.2828 + }
1.2829 +
1.2830 + /* === Generate the permutation ========================================= */
1.2831 +
1.2832 + for (c = 0 ; c < n_col ; c++)
1.2833 + {
1.2834 + p [Col [c].shared2.order] = c ;
1.2835 + }
1.2836 +}
1.2837 +
1.2838 +
1.2839 +/* ========================================================================== */
1.2840 +/* === detect_super_cols ==================================================== */
1.2841 +/* ========================================================================== */
1.2842 +
1.2843 +/*
1.2844 + Detects supercolumns by finding matches between columns in the hash buckets.
1.2845 + Check amongst columns in the set A [row_start ... row_start + row_length-1].
1.2846 + The columns under consideration are currently *not* in the degree lists,
1.2847 + and have already been placed in the hash buckets.
1.2848 +
1.2849 + The hash bucket for columns whose hash function is equal to h is stored
1.2850 + as follows:
1.2851 +
1.2852 + if head [h] is >= 0, then head [h] contains a degree list, so:
1.2853 +
1.2854 + head [h] is the first column in degree bucket h.
1.2855 + Col [head [h]].headhash gives the first column in hash bucket h.
1.2856 +
1.2857 + otherwise, the degree list is empty, and:
1.2858 +
1.2859 + -(head [h] + 2) is the first column in hash bucket h.
1.2860 +
1.2861 + For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
1.2862 + column" pointer. Col [c].shared3.hash is used instead as the hash number
1.2863 + for that column. The value of Col [c].shared4.hash_next is the next column
1.2864 + in the same hash bucket.
1.2865 +
1.2866 + Assuming no, or "few" hash collisions, the time taken by this routine is
1.2867 + linear in the sum of the sizes (lengths) of each column whose score has
1.2868 + just been computed in the approximate degree computation.
1.2869 + Not user-callable.
1.2870 +*/
1.2871 +
1.2872 +PRIVATE void detect_super_cols
1.2873 +(
1.2874 + /* === Parameters ======================================================= */
1.2875 +
1.2876 +#ifndef NDEBUG
1.2877 + /* these two parameters are only needed when debugging is enabled: */
1.2878 + Int n_col, /* number of columns of A */
1.2879 + Colamd_Row Row [], /* of size n_row+1 */
1.2880 +#endif /* NDEBUG */
1.2881 +
1.2882 + Colamd_Col Col [], /* of size n_col+1 */
1.2883 + Int A [], /* row indices of A */
1.2884 + Int head [], /* head of degree lists and hash buckets */
1.2885 + Int row_start, /* pointer to set of columns to check */
1.2886 + Int row_length /* number of columns to check */
1.2887 +)
1.2888 +{
1.2889 + /* === Local variables ================================================== */
1.2890 +
1.2891 + Int hash ; /* hash value for a column */
1.2892 + Int *rp ; /* pointer to a row */
1.2893 + Int c ; /* a column index */
1.2894 + Int super_c ; /* column index of the column to absorb into */
1.2895 + Int *cp1 ; /* column pointer for column super_c */
1.2896 + Int *cp2 ; /* column pointer for column c */
1.2897 + Int length ; /* length of column super_c */
1.2898 + Int prev_c ; /* column preceding c in hash bucket */
1.2899 + Int i ; /* loop counter */
1.2900 + Int *rp_end ; /* pointer to the end of the row */
1.2901 + Int col ; /* a column index in the row to check */
1.2902 + Int head_column ; /* first column in hash bucket or degree list */
1.2903 + Int first_col ; /* first column in hash bucket */
1.2904 +
1.2905 + /* === Consider each column in the row ================================== */
1.2906 +
1.2907 + rp = &A [row_start] ;
1.2908 + rp_end = rp + row_length ;
1.2909 + while (rp < rp_end)
1.2910 + {
1.2911 + col = *rp++ ;
1.2912 + if (COL_IS_DEAD (col))
1.2913 + {
1.2914 + continue ;
1.2915 + }
1.2916 +
1.2917 + /* get hash number for this column */
1.2918 + hash = Col [col].shared3.hash ;
1.2919 + ASSERT (hash <= n_col) ;
1.2920 +
1.2921 + /* === Get the first column in this hash bucket ===================== */
1.2922 +
1.2923 + head_column = head [hash] ;
1.2924 + if (head_column > EMPTY)
1.2925 + {
1.2926 + first_col = Col [head_column].shared3.headhash ;
1.2927 + }
1.2928 + else
1.2929 + {
1.2930 + first_col = - (head_column + 2) ;
1.2931 + }
1.2932 +
1.2933 + /* === Consider each column in the hash bucket ====================== */
1.2934 +
1.2935 + for (super_c = first_col ; super_c != EMPTY ;
1.2936 + super_c = Col [super_c].shared4.hash_next)
1.2937 + {
1.2938 + ASSERT (COL_IS_ALIVE (super_c)) ;
1.2939 + ASSERT (Col [super_c].shared3.hash == hash) ;
1.2940 + length = Col [super_c].length ;
1.2941 +
1.2942 + /* prev_c is the column preceding column c in the hash bucket */
1.2943 + prev_c = super_c ;
1.2944 +
1.2945 + /* === Compare super_c with all columns after it ================ */
1.2946 +
1.2947 + for (c = Col [super_c].shared4.hash_next ;
1.2948 + c != EMPTY ; c = Col [c].shared4.hash_next)
1.2949 + {
1.2950 + ASSERT (c != super_c) ;
1.2951 + ASSERT (COL_IS_ALIVE (c)) ;
1.2952 + ASSERT (Col [c].shared3.hash == hash) ;
1.2953 +
1.2954 + /* not identical if lengths or scores are different */
1.2955 + if (Col [c].length != length ||
1.2956 + Col [c].shared2.score != Col [super_c].shared2.score)
1.2957 + {
1.2958 + prev_c = c ;
1.2959 + continue ;
1.2960 + }
1.2961 +
1.2962 + /* compare the two columns */
1.2963 + cp1 = &A [Col [super_c].start] ;
1.2964 + cp2 = &A [Col [c].start] ;
1.2965 +
1.2966 + for (i = 0 ; i < length ; i++)
1.2967 + {
1.2968 + /* the columns are "clean" (no dead rows) */
1.2969 + ASSERT (ROW_IS_ALIVE (*cp1)) ;
1.2970 + ASSERT (ROW_IS_ALIVE (*cp2)) ;
1.2971 + /* row indices will same order for both supercols, */
1.2972 + /* no gather scatter nessasary */
1.2973 + if (*cp1++ != *cp2++)
1.2974 + {
1.2975 + break ;
1.2976 + }
1.2977 + }
1.2978 +
1.2979 + /* the two columns are different if the for-loop "broke" */
1.2980 + if (i != length)
1.2981 + {
1.2982 + prev_c = c ;
1.2983 + continue ;
1.2984 + }
1.2985 +
1.2986 + /* === Got it! two columns are identical =================== */
1.2987 +
1.2988 + ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
1.2989 +
1.2990 + Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
1.2991 + Col [c].shared1.parent = super_c ;
1.2992 + KILL_NON_PRINCIPAL_COL (c) ;
1.2993 + /* order c later, in order_children() */
1.2994 + Col [c].shared2.order = EMPTY ;
1.2995 + /* remove c from hash bucket */
1.2996 + Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
1.2997 + }
1.2998 + }
1.2999 +
1.3000 + /* === Empty this hash bucket ======================================= */
1.3001 +
1.3002 + if (head_column > EMPTY)
1.3003 + {
1.3004 + /* corresponding degree list "hash" is not empty */
1.3005 + Col [head_column].shared3.headhash = EMPTY ;
1.3006 + }
1.3007 + else
1.3008 + {
1.3009 + /* corresponding degree list "hash" is empty */
1.3010 + head [hash] = EMPTY ;
1.3011 + }
1.3012 + }
1.3013 +}
1.3014 +
1.3015 +
1.3016 +/* ========================================================================== */
1.3017 +/* === garbage_collection =================================================== */
1.3018 +/* ========================================================================== */
1.3019 +
1.3020 +/*
1.3021 + Defragments and compacts columns and rows in the workspace A. Used when
1.3022 + all avaliable memory has been used while performing row merging. Returns
1.3023 + the index of the first free position in A, after garbage collection. The
1.3024 + time taken by this routine is linear is the size of the array A, which is
1.3025 + itself linear in the number of nonzeros in the input matrix.
1.3026 + Not user-callable.
1.3027 +*/
1.3028 +
1.3029 +PRIVATE Int garbage_collection /* returns the new value of pfree */
1.3030 +(
1.3031 + /* === Parameters ======================================================= */
1.3032 +
1.3033 + Int n_row, /* number of rows */
1.3034 + Int n_col, /* number of columns */
1.3035 + Colamd_Row Row [], /* row info */
1.3036 + Colamd_Col Col [], /* column info */
1.3037 + Int A [], /* A [0 ... Alen-1] holds the matrix */
1.3038 + Int *pfree /* &A [0] ... pfree is in use */
1.3039 +)
1.3040 +{
1.3041 + /* === Local variables ================================================== */
1.3042 +
1.3043 + Int *psrc ; /* source pointer */
1.3044 + Int *pdest ; /* destination pointer */
1.3045 + Int j ; /* counter */
1.3046 + Int r ; /* a row index */
1.3047 + Int c ; /* a column index */
1.3048 + Int length ; /* length of a row or column */
1.3049 +
1.3050 +#ifndef NDEBUG
1.3051 + Int debug_rows ;
1.3052 + DEBUG2 (("Defrag..\n")) ;
1.3053 + for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
1.3054 + debug_rows = 0 ;
1.3055 +#endif /* NDEBUG */
1.3056 +
1.3057 + /* === Defragment the columns =========================================== */
1.3058 +
1.3059 + pdest = &A[0] ;
1.3060 + for (c = 0 ; c < n_col ; c++)
1.3061 + {
1.3062 + if (COL_IS_ALIVE (c))
1.3063 + {
1.3064 + psrc = &A [Col [c].start] ;
1.3065 +
1.3066 + /* move and compact the column */
1.3067 + ASSERT (pdest <= psrc) ;
1.3068 + Col [c].start = (Int) (pdest - &A [0]) ;
1.3069 + length = Col [c].length ;
1.3070 + for (j = 0 ; j < length ; j++)
1.3071 + {
1.3072 + r = *psrc++ ;
1.3073 + if (ROW_IS_ALIVE (r))
1.3074 + {
1.3075 + *pdest++ = r ;
1.3076 + }
1.3077 + }
1.3078 + Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
1.3079 + }
1.3080 + }
1.3081 +
1.3082 + /* === Prepare to defragment the rows =================================== */
1.3083 +
1.3084 + for (r = 0 ; r < n_row ; r++)
1.3085 + {
1.3086 + if (ROW_IS_DEAD (r) || (Row [r].length == 0))
1.3087 + {
1.3088 + /* This row is already dead, or is of zero length. Cannot compact
1.3089 + * a row of zero length, so kill it. NOTE: in the current version,
1.3090 + * there are no zero-length live rows. Kill the row (for the first
1.3091 + * time, or again) just to be safe. */
1.3092 + KILL_ROW (r) ;
1.3093 + }
1.3094 + else
1.3095 + {
1.3096 + /* save first column index in Row [r].shared2.first_column */
1.3097 + psrc = &A [Row [r].start] ;
1.3098 + Row [r].shared2.first_column = *psrc ;
1.3099 + ASSERT (ROW_IS_ALIVE (r)) ;
1.3100 + /* flag the start of the row with the one's complement of row */
1.3101 + *psrc = ONES_COMPLEMENT (r) ;
1.3102 +#ifndef NDEBUG
1.3103 + debug_rows++ ;
1.3104 +#endif /* NDEBUG */
1.3105 + }
1.3106 + }
1.3107 +
1.3108 + /* === Defragment the rows ============================================== */
1.3109 +
1.3110 + psrc = pdest ;
1.3111 + while (psrc < pfree)
1.3112 + {
1.3113 + /* find a negative number ... the start of a row */
1.3114 + if (*psrc++ < 0)
1.3115 + {
1.3116 + psrc-- ;
1.3117 + /* get the row index */
1.3118 + r = ONES_COMPLEMENT (*psrc) ;
1.3119 + ASSERT (r >= 0 && r < n_row) ;
1.3120 + /* restore first column index */
1.3121 + *psrc = Row [r].shared2.first_column ;
1.3122 + ASSERT (ROW_IS_ALIVE (r)) ;
1.3123 + ASSERT (Row [r].length > 0) ;
1.3124 + /* move and compact the row */
1.3125 + ASSERT (pdest <= psrc) ;
1.3126 + Row [r].start = (Int) (pdest - &A [0]) ;
1.3127 + length = Row [r].length ;
1.3128 + for (j = 0 ; j < length ; j++)
1.3129 + {
1.3130 + c = *psrc++ ;
1.3131 + if (COL_IS_ALIVE (c))
1.3132 + {
1.3133 + *pdest++ = c ;
1.3134 + }
1.3135 + }
1.3136 + Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
1.3137 + ASSERT (Row [r].length > 0) ;
1.3138 +#ifndef NDEBUG
1.3139 + debug_rows-- ;
1.3140 +#endif /* NDEBUG */
1.3141 + }
1.3142 + }
1.3143 + /* ensure we found all the rows */
1.3144 + ASSERT (debug_rows == 0) ;
1.3145 +
1.3146 + /* === Return the new value of pfree ==================================== */
1.3147 +
1.3148 + return ((Int) (pdest - &A [0])) ;
1.3149 +}
1.3150 +
1.3151 +
1.3152 +/* ========================================================================== */
1.3153 +/* === clear_mark =========================================================== */
1.3154 +/* ========================================================================== */
1.3155 +
1.3156 +/*
1.3157 + Clears the Row [].shared2.mark array, and returns the new tag_mark.
1.3158 + Return value is the new tag_mark. Not user-callable.
1.3159 +*/
1.3160 +
1.3161 +PRIVATE Int clear_mark /* return the new value for tag_mark */
1.3162 +(
1.3163 + /* === Parameters ======================================================= */
1.3164 +
1.3165 + Int tag_mark, /* new value of tag_mark */
1.3166 + Int max_mark, /* max allowed value of tag_mark */
1.3167 +
1.3168 + Int n_row, /* number of rows in A */
1.3169 + Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
1.3170 +)
1.3171 +{
1.3172 + /* === Local variables ================================================== */
1.3173 +
1.3174 + Int r ;
1.3175 +
1.3176 + if (tag_mark <= 0 || tag_mark >= max_mark)
1.3177 + {
1.3178 + for (r = 0 ; r < n_row ; r++)
1.3179 + {
1.3180 + if (ROW_IS_ALIVE (r))
1.3181 + {
1.3182 + Row [r].shared2.mark = 0 ;
1.3183 + }
1.3184 + }
1.3185 + tag_mark = 1 ;
1.3186 + }
1.3187 +
1.3188 + return (tag_mark) ;
1.3189 +}
1.3190 +
1.3191 +
1.3192 +/* ========================================================================== */
1.3193 +/* === print_report ========================================================= */
1.3194 +/* ========================================================================== */
1.3195 +
1.3196 +PRIVATE void print_report
1.3197 +(
1.3198 + char *method,
1.3199 + Int stats [COLAMD_STATS]
1.3200 +)
1.3201 +{
1.3202 +
1.3203 + Int i1, i2, i3 ;
1.3204 +
1.3205 + PRINTF (("\n%s version %d.%d, %s: ", method,
1.3206 + COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
1.3207 +
1.3208 + if (!stats)
1.3209 + {
1.3210 + PRINTF (("No statistics available.\n")) ;
1.3211 + return ;
1.3212 + }
1.3213 +
1.3214 + i1 = stats [COLAMD_INFO1] ;
1.3215 + i2 = stats [COLAMD_INFO2] ;
1.3216 + i3 = stats [COLAMD_INFO3] ;
1.3217 +
1.3218 + if (stats [COLAMD_STATUS] >= 0)
1.3219 + {
1.3220 + PRINTF (("OK. ")) ;
1.3221 + }
1.3222 + else
1.3223 + {
1.3224 + PRINTF (("ERROR. ")) ;
1.3225 + }
1.3226 +
1.3227 + switch (stats [COLAMD_STATUS])
1.3228 + {
1.3229 +
1.3230 + case COLAMD_OK_BUT_JUMBLED:
1.3231 +
1.3232 + PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
1.3233 +
1.3234 + PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
1.3235 + method, i3)) ;
1.3236 +
1.3237 + PRINTF(("%s: last seen duplicate or out-of-order row index: %d\n",
1.3238 + method, INDEX (i2))) ;
1.3239 +
1.3240 + PRINTF(("%s: last seen in column: %d",
1.3241 + method, INDEX (i1))) ;
1.3242 +
1.3243 + /* no break - fall through to next case instead */
1.3244 +
1.3245 + case COLAMD_OK:
1.3246 +
1.3247 + PRINTF(("\n")) ;
1.3248 +
1.3249 + PRINTF(("%s: number of dense or empty rows ignored: %d\n",
1.3250 + method, stats [COLAMD_DENSE_ROW])) ;
1.3251 +
1.3252 + PRINTF(("%s: number of dense or empty columns ignored: %d\n",
1.3253 + method, stats [COLAMD_DENSE_COL])) ;
1.3254 +
1.3255 + PRINTF(("%s: number of garbage collections performed: %d\n",
1.3256 + method, stats [COLAMD_DEFRAG_COUNT])) ;
1.3257 + break ;
1.3258 +
1.3259 + case COLAMD_ERROR_A_not_present:
1.3260 +
1.3261 + PRINTF(("Array A (row indices of matrix) not present.\n")) ;
1.3262 + break ;
1.3263 +
1.3264 + case COLAMD_ERROR_p_not_present:
1.3265 +
1.3266 + PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
1.3267 + break ;
1.3268 +
1.3269 + case COLAMD_ERROR_nrow_negative:
1.3270 +
1.3271 + PRINTF(("Invalid number of rows (%d).\n", i1)) ;
1.3272 + break ;
1.3273 +
1.3274 + case COLAMD_ERROR_ncol_negative:
1.3275 +
1.3276 + PRINTF(("Invalid number of columns (%d).\n", i1)) ;
1.3277 + break ;
1.3278 +
1.3279 + case COLAMD_ERROR_nnz_negative:
1.3280 +
1.3281 + PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
1.3282 + break ;
1.3283 +
1.3284 + case COLAMD_ERROR_p0_nonzero:
1.3285 +
1.3286 + PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
1.3287 + break ;
1.3288 +
1.3289 + case COLAMD_ERROR_A_too_small:
1.3290 +
1.3291 + PRINTF(("Array A too small.\n")) ;
1.3292 + PRINTF((" Need Alen >= %d, but given only Alen = %d.\n",
1.3293 + i1, i2)) ;
1.3294 + break ;
1.3295 +
1.3296 + case COLAMD_ERROR_col_length_negative:
1.3297 +
1.3298 + PRINTF
1.3299 + (("Column %d has a negative number of nonzero entries (%d).\n",
1.3300 + INDEX (i1), i2)) ;
1.3301 + break ;
1.3302 +
1.3303 + case COLAMD_ERROR_row_index_out_of_bounds:
1.3304 +
1.3305 + PRINTF
1.3306 + (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
1.3307 + INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
1.3308 + break ;
1.3309 +
1.3310 + case COLAMD_ERROR_out_of_memory:
1.3311 +
1.3312 + PRINTF(("Out of memory.\n")) ;
1.3313 + break ;
1.3314 +
1.3315 + /* v2.4: internal-error case deleted */
1.3316 + }
1.3317 +}
1.3318 +
1.3319 +
1.3320 +
1.3321 +
1.3322 +/* ========================================================================== */
1.3323 +/* === colamd debugging routines ============================================ */
1.3324 +/* ========================================================================== */
1.3325 +
1.3326 +/* When debugging is disabled, the remainder of this file is ignored. */
1.3327 +
1.3328 +#ifndef NDEBUG
1.3329 +
1.3330 +
1.3331 +/* ========================================================================== */
1.3332 +/* === debug_structures ===================================================== */
1.3333 +/* ========================================================================== */
1.3334 +
1.3335 +/*
1.3336 + At this point, all empty rows and columns are dead. All live columns
1.3337 + are "clean" (containing no dead rows) and simplicial (no supercolumns
1.3338 + yet). Rows may contain dead columns, but all live rows contain at
1.3339 + least one live column.
1.3340 +*/
1.3341 +
1.3342 +PRIVATE void debug_structures
1.3343 +(
1.3344 + /* === Parameters ======================================================= */
1.3345 +
1.3346 + Int n_row,
1.3347 + Int n_col,
1.3348 + Colamd_Row Row [],
1.3349 + Colamd_Col Col [],
1.3350 + Int A [],
1.3351 + Int n_col2
1.3352 +)
1.3353 +{
1.3354 + /* === Local variables ================================================== */
1.3355 +
1.3356 + Int i ;
1.3357 + Int c ;
1.3358 + Int *cp ;
1.3359 + Int *cp_end ;
1.3360 + Int len ;
1.3361 + Int score ;
1.3362 + Int r ;
1.3363 + Int *rp ;
1.3364 + Int *rp_end ;
1.3365 + Int deg ;
1.3366 +
1.3367 + /* === Check A, Row, and Col ============================================ */
1.3368 +
1.3369 + for (c = 0 ; c < n_col ; c++)
1.3370 + {
1.3371 + if (COL_IS_ALIVE (c))
1.3372 + {
1.3373 + len = Col [c].length ;
1.3374 + score = Col [c].shared2.score ;
1.3375 + DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
1.3376 + ASSERT (len > 0) ;
1.3377 + ASSERT (score >= 0) ;
1.3378 + ASSERT (Col [c].shared1.thickness == 1) ;
1.3379 + cp = &A [Col [c].start] ;
1.3380 + cp_end = cp + len ;
1.3381 + while (cp < cp_end)
1.3382 + {
1.3383 + r = *cp++ ;
1.3384 + ASSERT (ROW_IS_ALIVE (r)) ;
1.3385 + }
1.3386 + }
1.3387 + else
1.3388 + {
1.3389 + i = Col [c].shared2.order ;
1.3390 + ASSERT (i >= n_col2 && i < n_col) ;
1.3391 + }
1.3392 + }
1.3393 +
1.3394 + for (r = 0 ; r < n_row ; r++)
1.3395 + {
1.3396 + if (ROW_IS_ALIVE (r))
1.3397 + {
1.3398 + i = 0 ;
1.3399 + len = Row [r].length ;
1.3400 + deg = Row [r].shared1.degree ;
1.3401 + ASSERT (len > 0) ;
1.3402 + ASSERT (deg > 0) ;
1.3403 + rp = &A [Row [r].start] ;
1.3404 + rp_end = rp + len ;
1.3405 + while (rp < rp_end)
1.3406 + {
1.3407 + c = *rp++ ;
1.3408 + if (COL_IS_ALIVE (c))
1.3409 + {
1.3410 + i++ ;
1.3411 + }
1.3412 + }
1.3413 + ASSERT (i > 0) ;
1.3414 + }
1.3415 + }
1.3416 +}
1.3417 +
1.3418 +
1.3419 +/* ========================================================================== */
1.3420 +/* === debug_deg_lists ====================================================== */
1.3421 +/* ========================================================================== */
1.3422 +
1.3423 +/*
1.3424 + Prints the contents of the degree lists. Counts the number of columns
1.3425 + in the degree list and compares it to the total it should have. Also
1.3426 + checks the row degrees.
1.3427 +*/
1.3428 +
1.3429 +PRIVATE void debug_deg_lists
1.3430 +(
1.3431 + /* === Parameters ======================================================= */
1.3432 +
1.3433 + Int n_row,
1.3434 + Int n_col,
1.3435 + Colamd_Row Row [],
1.3436 + Colamd_Col Col [],
1.3437 + Int head [],
1.3438 + Int min_score,
1.3439 + Int should,
1.3440 + Int max_deg
1.3441 +)
1.3442 +{
1.3443 + /* === Local variables ================================================== */
1.3444 +
1.3445 + Int deg ;
1.3446 + Int col ;
1.3447 + Int have ;
1.3448 + Int row ;
1.3449 +
1.3450 + /* === Check the degree lists =========================================== */
1.3451 +
1.3452 + if (n_col > 10000 && colamd_debug <= 0)
1.3453 + {
1.3454 + return ;
1.3455 + }
1.3456 + have = 0 ;
1.3457 + DEBUG4 (("Degree lists: %d\n", min_score)) ;
1.3458 + for (deg = 0 ; deg <= n_col ; deg++)
1.3459 + {
1.3460 + col = head [deg] ;
1.3461 + if (col == EMPTY)
1.3462 + {
1.3463 + continue ;
1.3464 + }
1.3465 + DEBUG4 (("%d:", deg)) ;
1.3466 + while (col != EMPTY)
1.3467 + {
1.3468 + DEBUG4 ((" %d", col)) ;
1.3469 + have += Col [col].shared1.thickness ;
1.3470 + ASSERT (COL_IS_ALIVE (col)) ;
1.3471 + col = Col [col].shared4.degree_next ;
1.3472 + }
1.3473 + DEBUG4 (("\n")) ;
1.3474 + }
1.3475 + DEBUG4 (("should %d have %d\n", should, have)) ;
1.3476 + ASSERT (should == have) ;
1.3477 +
1.3478 + /* === Check the row degrees ============================================ */
1.3479 +
1.3480 + if (n_row > 10000 && colamd_debug <= 0)
1.3481 + {
1.3482 + return ;
1.3483 + }
1.3484 + for (row = 0 ; row < n_row ; row++)
1.3485 + {
1.3486 + if (ROW_IS_ALIVE (row))
1.3487 + {
1.3488 + ASSERT (Row [row].shared1.degree <= max_deg) ;
1.3489 + }
1.3490 + }
1.3491 +}
1.3492 +
1.3493 +
1.3494 +/* ========================================================================== */
1.3495 +/* === debug_mark =========================================================== */
1.3496 +/* ========================================================================== */
1.3497 +
1.3498 +/*
1.3499 + Ensures that the tag_mark is less that the maximum and also ensures that
1.3500 + each entry in the mark array is less than the tag mark.
1.3501 +*/
1.3502 +
1.3503 +PRIVATE void debug_mark
1.3504 +(
1.3505 + /* === Parameters ======================================================= */
1.3506 +
1.3507 + Int n_row,
1.3508 + Colamd_Row Row [],
1.3509 + Int tag_mark,
1.3510 + Int max_mark
1.3511 +)
1.3512 +{
1.3513 + /* === Local variables ================================================== */
1.3514 +
1.3515 + Int r ;
1.3516 +
1.3517 + /* === Check the Row marks ============================================== */
1.3518 +
1.3519 + ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
1.3520 + if (n_row > 10000 && colamd_debug <= 0)
1.3521 + {
1.3522 + return ;
1.3523 + }
1.3524 + for (r = 0 ; r < n_row ; r++)
1.3525 + {
1.3526 + ASSERT (Row [r].shared2.mark < tag_mark) ;
1.3527 + }
1.3528 +}
1.3529 +
1.3530 +
1.3531 +/* ========================================================================== */
1.3532 +/* === debug_matrix ========================================================= */
1.3533 +/* ========================================================================== */
1.3534 +
1.3535 +/*
1.3536 + Prints out the contents of the columns and the rows.
1.3537 +*/
1.3538 +
1.3539 +PRIVATE void debug_matrix
1.3540 +(
1.3541 + /* === Parameters ======================================================= */
1.3542 +
1.3543 + Int n_row,
1.3544 + Int n_col,
1.3545 + Colamd_Row Row [],
1.3546 + Colamd_Col Col [],
1.3547 + Int A []
1.3548 +)
1.3549 +{
1.3550 + /* === Local variables ================================================== */
1.3551 +
1.3552 + Int r ;
1.3553 + Int c ;
1.3554 + Int *rp ;
1.3555 + Int *rp_end ;
1.3556 + Int *cp ;
1.3557 + Int *cp_end ;
1.3558 +
1.3559 + /* === Dump the rows and columns of the matrix ========================== */
1.3560 +
1.3561 + if (colamd_debug < 3)
1.3562 + {
1.3563 + return ;
1.3564 + }
1.3565 + DEBUG3 (("DUMP MATRIX:\n")) ;
1.3566 + for (r = 0 ; r < n_row ; r++)
1.3567 + {
1.3568 + DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
1.3569 + if (ROW_IS_DEAD (r))
1.3570 + {
1.3571 + continue ;
1.3572 + }
1.3573 + DEBUG3 (("start %d length %d degree %d\n",
1.3574 + Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
1.3575 + rp = &A [Row [r].start] ;
1.3576 + rp_end = rp + Row [r].length ;
1.3577 + while (rp < rp_end)
1.3578 + {
1.3579 + c = *rp++ ;
1.3580 + DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ;
1.3581 + }
1.3582 + }
1.3583 +
1.3584 + for (c = 0 ; c < n_col ; c++)
1.3585 + {
1.3586 + DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
1.3587 + if (COL_IS_DEAD (c))
1.3588 + {
1.3589 + continue ;
1.3590 + }
1.3591 + DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
1.3592 + Col [c].start, Col [c].length,
1.3593 + Col [c].shared1.thickness, Col [c].shared2.score)) ;
1.3594 + cp = &A [Col [c].start] ;
1.3595 + cp_end = cp + Col [c].length ;
1.3596 + while (cp < cp_end)
1.3597 + {
1.3598 + r = *cp++ ;
1.3599 + DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ;
1.3600 + }
1.3601 + }
1.3602 +}
1.3603 +
1.3604 +PRIVATE void colamd_get_debug
1.3605 +(
1.3606 + char *method
1.3607 +)
1.3608 +{
1.3609 + FILE *f ;
1.3610 + colamd_debug = 0 ; /* no debug printing */
1.3611 + f = fopen ("debug", "r") ;
1.3612 + if (f == (FILE *) NULL)
1.3613 + {
1.3614 + colamd_debug = 0 ;
1.3615 + }
1.3616 + else
1.3617 + {
1.3618 + fscanf (f, "%d", &colamd_debug) ;
1.3619 + fclose (f) ;
1.3620 + }
1.3621 + DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
1.3622 + method, colamd_debug)) ;
1.3623 +}
1.3624 +
1.3625 +#endif /* NDEBUG */