lemon-project-template-glpk: deps/glpk/src/colamd/colamd.c annotate

lemon-project-template-glpk

annotate deps/glpk/src/colamd/colamd.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

rev	line source
alpar@9	1 /* ========================================================================== */
alpar@9	2 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
alpar@9	3 /* ========================================================================== */
alpar@9	4
alpar@9	5 /* COLAMD / SYMAMD
alpar@9	6
alpar@9	7 colamd: an approximate minimum degree column ordering algorithm,
alpar@9	8 for LU factorization of symmetric or unsymmetric matrices,
alpar@9	9 QR factorization, least squares, interior point methods for
alpar@9	10 linear programming problems, and other related problems.
alpar@9	11
alpar@9	12 symamd: an approximate minimum degree ordering algorithm for Cholesky
alpar@9	13 factorization of symmetric matrices.
alpar@9	14
alpar@9	15 Purpose:
alpar@9	16
alpar@9	17 Colamd computes a permutation Q such that the Cholesky factorization of
alpar@9	18 (AQ)'(AQ) has less fill-in and requires fewer floating point operations
alpar@9	19 than A'A. This also provides a good ordering for sparse partial
alpar@9	20 pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
alpar@9	21 factorization, and P is computed during numerical factorization via
alpar@9	22 conventional partial pivoting with row interchanges. Colamd is the
alpar@9	23 column ordering method used in SuperLU, part of the ScaLAPACK library.
alpar@9	24 It is also available as built-in function in MATLAB Version 6,
alpar@9	25 available from MathWorks, Inc. (http://www.mathworks.com). This
alpar@9	26 routine can be used in place of colmmd in MATLAB.
alpar@9	27
alpar@9	28 Symamd computes a permutation P of a symmetric matrix A such that the
alpar@9	29 Cholesky factorization of PAP' has less fill-in and requires fewer
alpar@9	30 floating point operations than A. Symamd constructs a matrix M such
alpar@9	31 that M'M has the same nonzero pattern of A, and then orders the columns
alpar@9	32 of M using colmmd. The column ordering of M is then returned as the
alpar@9	33 row and column ordering P of A.
alpar@9	34
alpar@9	35 Authors:
alpar@9	36
alpar@9	37 The authors of the code itself are Stefan I. Larimore and Timothy A.
alpar@9	38 Davis (davis at cise.ufl.edu), University of Florida. The algorithm was
alpar@9	39 developed in collaboration with John Gilbert, Xerox PARC, and Esmond
alpar@9	40 Ng, Oak Ridge National Laboratory.
alpar@9	41
alpar@9	42 Acknowledgements:
alpar@9	43
alpar@9	44 This work was supported by the National Science Foundation, under
alpar@9	45 grants DMS-9504974 and DMS-9803599.
alpar@9	46
alpar@9	47 Copyright and License:
alpar@9	48
alpar@9	49 Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
alpar@9	50 COLAMD is also available under alternate licenses, contact T. Davis
alpar@9	51 for details.
alpar@9	52
alpar@9	53 This library is free software; you can redistribute it and/or
alpar@9	54 modify it under the terms of the GNU Lesser General Public
alpar@9	55 License as published by the Free Software Foundation; either
alpar@9	56 version 2.1 of the License, or (at your option) any later version.
alpar@9	57
alpar@9	58 This library is distributed in the hope that it will be useful,
alpar@9	59 but WITHOUT ANY WARRANTY; without even the implied warranty of
alpar@9	60 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
alpar@9	61 Lesser General Public License for more details.
alpar@9	62
alpar@9	63 You should have received a copy of the GNU Lesser General Public
alpar@9	64 License along with this library; if not, write to the Free Software
alpar@9	65 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
alpar@9	66 USA
alpar@9	67
alpar@9	68 Permission is hereby granted to use or copy this program under the
alpar@9	69 terms of the GNU LGPL, provided that the Copyright, this License,
alpar@9	70 and the Availability of the original version is retained on all copies.
alpar@9	71 User documentation of any code that uses this code or any modified
alpar@9	72 version of this code must cite the Copyright, this License, the
alpar@9	73 Availability note, and "Used by permission." Permission to modify
alpar@9	74 the code and to distribute modified code is granted, provided the
alpar@9	75 Copyright, this License, and the Availability note are retained,
alpar@9	76 and a notice that the code was modified is included.
alpar@9	77
alpar@9	78 Availability:
alpar@9	79
alpar@9	80 The colamd/symamd library is available at
alpar@9	81
alpar@9	82 http://www.cise.ufl.edu/research/sparse/colamd/
alpar@9	83
alpar@9	84 This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
alpar@9	85 file. It requires the colamd.h file. It is required by the colamdmex.c
alpar@9	86 and symamdmex.c files, for the MATLAB interface to colamd and symamd.
alpar@9	87 Appears as ACM Algorithm 836.
alpar@9	88
alpar@9	89 See the ChangeLog file for changes since Version 1.0.
alpar@9	90
alpar@9	91 References:
alpar@9	92
alpar@9	93 T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
alpar@9	94 minimum degree ordering algorithm, ACM Transactions on Mathematical
alpar@9	95 Software, vol. 30, no. 3., pp. 353-376, 2004.
alpar@9	96
alpar@9	97 T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
alpar@9	98 an approximate column minimum degree ordering algorithm, ACM
alpar@9	99 Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
alpar@9	100 2004.
alpar@9	101
alpar@9	102 */
alpar@9	103
alpar@9	104 /* ========================================================================== */
alpar@9	105 /* === Description of user-callable routines ================================ */
alpar@9	106 /* ========================================================================== */
alpar@9	107
alpar@9	108 /* COLAMD includes both int and UF_long versions of all its routines. The
alpar@9	109 * description below is for the int version. For UF_long, all int arguments
alpar@9	110 * become UF_long. UF_long is normally defined as long, except for WIN64.
alpar@9	111
alpar@9	112 ----------------------------------------------------------------------------
alpar@9	113 colamd_recommended:
alpar@9	114 ----------------------------------------------------------------------------
alpar@9	115
alpar@9	116 C syntax:
alpar@9	117
alpar@9	118 #include "colamd.h"
alpar@9	119 size_t colamd_recommended (int nnz, int n_row, int n_col) ;
alpar@9	120 size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
alpar@9	121 UF_long n_col) ;
alpar@9	122
alpar@9	123 Purpose:
alpar@9	124
alpar@9	125 Returns recommended value of Alen for use by colamd. Returns 0
alpar@9	126 if any input argument is negative. The use of this routine
alpar@9	127 is optional. Not needed for symamd, which dynamically allocates
alpar@9	128 its own memory.
alpar@9	129
alpar@9	130 Note that in v2.4 and earlier, these routines returned int or long.
alpar@9	131 They now return a value of type size_t.
alpar@9	132
alpar@9	133 Arguments (all input arguments):
alpar@9	134
alpar@9	135 int nnz ; Number of nonzeros in the matrix A. This must
alpar@9	136 be the same value as p [n_col] in the call to
alpar@9	137 colamd - otherwise you will get a wrong value
alpar@9	138 of the recommended memory to use.
alpar@9	139
alpar@9	140 int n_row ; Number of rows in the matrix A.
alpar@9	141
alpar@9	142 int n_col ; Number of columns in the matrix A.
alpar@9	143
alpar@9	144 ----------------------------------------------------------------------------
alpar@9	145 colamd_set_defaults:
alpar@9	146 ----------------------------------------------------------------------------
alpar@9	147
alpar@9	148 C syntax:
alpar@9	149
alpar@9	150 #include "colamd.h"
alpar@9	151 colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
alpar@9	152 colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
alpar@9	153
alpar@9	154 Purpose:
alpar@9	155
alpar@9	156 Sets the default parameters. The use of this routine is optional.
alpar@9	157
alpar@9	158 Arguments:
alpar@9	159
alpar@9	160 double knobs [COLAMD_KNOBS] ; Output only.
alpar@9	161
alpar@9	162 NOTE: the meaning of the dense row/col knobs has changed in v2.4
alpar@9	163
alpar@9	164 knobs [0] and knobs [1] control dense row and col detection:
alpar@9	165
alpar@9	166 Colamd: rows with more than
alpar@9	167 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
alpar@9	168 entries are removed prior to ordering. Columns with more than
alpar@9	169 max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
alpar@9	170 entries are removed prior to
alpar@9	171 ordering, and placed last in the output column ordering.
alpar@9	172
alpar@9	173 Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
alpar@9	174 Rows and columns with more than
alpar@9	175 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
alpar@9	176 entries are removed prior to ordering, and placed last in the
alpar@9	177 output ordering.
alpar@9	178
alpar@9	179 COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
alpar@9	180 respectively, in colamd.h. Default values of these two knobs
alpar@9	181 are both 10. Currently, only knobs [0] and knobs [1] are
alpar@9	182 used, but future versions may use more knobs. If so, they will
alpar@9	183 be properly set to their defaults by the future version of
alpar@9	184 colamd_set_defaults, so that the code that calls colamd will
alpar@9	185 not need to change, assuming that you either use
alpar@9	186 colamd_set_defaults, or pass a (double *) NULL pointer as the
alpar@9	187 knobs array to colamd or symamd.
alpar@9	188
alpar@9	189 knobs [2]: aggressive absorption
alpar@9	190
alpar@9	191 knobs [COLAMD_AGGRESSIVE] controls whether or not to do
alpar@9	192 aggressive absorption during the ordering. Default is TRUE.
alpar@9	193
alpar@9	194
alpar@9	195 ----------------------------------------------------------------------------
alpar@9	196 colamd:
alpar@9	197 ----------------------------------------------------------------------------
alpar@9	198
alpar@9	199 C syntax:
alpar@9	200
alpar@9	201 #include "colamd.h"
alpar@9	202 int colamd (int n_row, int n_col, int Alen, int A, int p,
alpar@9	203 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
alpar@9	204 UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
alpar@9	205 UF_long A, UF_long p, double knobs [COLAMD_KNOBS],
alpar@9	206 UF_long stats [COLAMD_STATS]) ;
alpar@9	207
alpar@9	208 Purpose:
alpar@9	209
alpar@9	210 Computes a column ordering (Q) of A such that P(AQ)=LU or
alpar@9	211 (AQ)'AQ=LL' have less fill-in and require fewer floating point
alpar@9	212 operations than factorizing the unpermuted matrix A or A'A,
alpar@9	213 respectively.
alpar@9	214
alpar@9	215 Returns:
alpar@9	216
alpar@9	217 TRUE (1) if successful, FALSE (0) otherwise.
alpar@9	218
alpar@9	219 Arguments:
alpar@9	220
alpar@9	221 int n_row ; Input argument.
alpar@9	222
alpar@9	223 Number of rows in the matrix A.
alpar@9	224 Restriction: n_row >= 0.
alpar@9	225 Colamd returns FALSE if n_row is negative.
alpar@9	226
alpar@9	227 int n_col ; Input argument.
alpar@9	228
alpar@9	229 Number of columns in the matrix A.
alpar@9	230 Restriction: n_col >= 0.
alpar@9	231 Colamd returns FALSE if n_col is negative.
alpar@9	232
alpar@9	233 int Alen ; Input argument.
alpar@9	234
alpar@9	235 Restriction (see note):
alpar@9	236 Alen >= 2nnz + 6(n_col+1) + 4*(n_row+1) + n_col
alpar@9	237 Colamd returns FALSE if these conditions are not met.
alpar@9	238
alpar@9	239 Note: this restriction makes an modest assumption regarding
alpar@9	240 the size of the two typedef's structures in colamd.h.
alpar@9	241 We do, however, guarantee that
alpar@9	242
alpar@9	243 Alen >= colamd_recommended (nnz, n_row, n_col)
alpar@9	244
alpar@9	245 will be sufficient. Note: the macro version does not check
alpar@9	246 for integer overflow, and thus is not recommended. Use
alpar@9	247 the colamd_recommended routine instead.
alpar@9	248
alpar@9	249 int A [Alen] ; Input argument, undefined on output.
alpar@9	250
alpar@9	251 A is an integer array of size Alen. Alen must be at least as
alpar@9	252 large as the bare minimum value given above, but this is very
alpar@9	253 low, and can result in excessive run time. For best
alpar@9	254 performance, we recommend that Alen be greater than or equal to
alpar@9	255 colamd_recommended (nnz, n_row, n_col), which adds
alpar@9	256 nnz/5 to the bare minimum value given above.
alpar@9	257
alpar@9	258 On input, the row indices of the entries in column c of the
alpar@9	259 matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices
alpar@9	260 in a given column c need not be in ascending order, and
alpar@9	261 duplicate row indices may be be present. However, colamd will
alpar@9	262 work a little faster if both of these conditions are met
alpar@9	263 (Colamd puts the matrix into this format, if it finds that the
alpar@9	264 the conditions are not met).
alpar@9	265
alpar@9	266 The matrix is 0-based. That is, rows are in the range 0 to
alpar@9	267 n_row-1, and columns are in the range 0 to n_col-1. Colamd
alpar@9	268 returns FALSE if any row index is out of range.
alpar@9	269
alpar@9	270 The contents of A are modified during ordering, and are
alpar@9	271 undefined on output.
alpar@9	272
alpar@9	273 int p [n_col+1] ; Both input and output argument.
alpar@9	274
alpar@9	275 p is an integer array of size n_col+1. On input, it holds the
alpar@9	276 "pointers" for the column form of the matrix A. Column c of
alpar@9	277 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
alpar@9	278 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
alpar@9	279 for all c in the range 0 to n_col-1. The value p [n_col] is
alpar@9	280 thus the total number of entries in the pattern of the matrix A.
alpar@9	281 Colamd returns FALSE if these conditions are not met.
alpar@9	282
alpar@9	283 On output, if colamd returns TRUE, the array p holds the column
alpar@9	284 permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
alpar@9	285 the first column index in the new ordering, and p [n_col-1] is
alpar@9	286 the last. That is, p [k] = j means that column j of A is the
alpar@9	287 kth pivot column, in AQ, where k is in the range 0 to n_col-1
alpar@9	288 (p [0] = j means that column j of A is the first column in AQ).
alpar@9	289
alpar@9	290 If colamd returns FALSE, then no permutation is returned, and
alpar@9	291 p is undefined on output.
alpar@9	292
alpar@9	293 double knobs [COLAMD_KNOBS] ; Input argument.
alpar@9	294
alpar@9	295 See colamd_set_defaults for a description.
alpar@9	296
alpar@9	297 int stats [COLAMD_STATS] ; Output argument.
alpar@9	298
alpar@9	299 Statistics on the ordering, and error status.
alpar@9	300 See colamd.h for related definitions.
alpar@9	301 Colamd returns FALSE if stats is not present.
alpar@9	302
alpar@9	303 stats [0]: number of dense or empty rows ignored.
alpar@9	304
alpar@9	305 stats [1]: number of dense or empty columns ignored (and
alpar@9	306 ordered last in the output permutation p)
alpar@9	307 Note that a row can become "empty" if it
alpar@9	308 contains only "dense" and/or "empty" columns,
alpar@9	309 and similarly a column can become "empty" if it
alpar@9	310 only contains "dense" and/or "empty" rows.
alpar@9	311
alpar@9	312 stats [2]: number of garbage collections performed.
alpar@9	313 This can be excessively high if Alen is close
alpar@9	314 to the minimum required value.
alpar@9	315
alpar@9	316 stats [3]: status code. < 0 is an error code.
alpar@9	317 > 1 is a warning or notice.
alpar@9	318
alpar@9	319 0 OK. Each column of the input matrix contained
alpar@9	320 row indices in increasing order, with no
alpar@9	321 duplicates.
alpar@9	322
alpar@9	323 1 OK, but columns of input matrix were jumbled
alpar@9	324 (unsorted columns or duplicate entries). Colamd
alpar@9	325 had to do some extra work to sort the matrix
alpar@9	326 first and remove duplicate entries, but it
alpar@9	327 still was able to return a valid permutation
alpar@9	328 (return value of colamd was TRUE).
alpar@9	329
alpar@9	330 stats [4]: highest numbered column that
alpar@9	331 is unsorted or has duplicate
alpar@9	332 entries.
alpar@9	333 stats [5]: last seen duplicate or
alpar@9	334 unsorted row index.
alpar@9	335 stats [6]: number of duplicate or
alpar@9	336 unsorted row indices.
alpar@9	337
alpar@9	338 -1 A is a null pointer
alpar@9	339
alpar@9	340 -2 p is a null pointer
alpar@9	341
alpar@9	342 -3 n_row is negative
alpar@9	343
alpar@9	344 stats [4]: n_row
alpar@9	345
alpar@9	346 -4 n_col is negative
alpar@9	347
alpar@9	348 stats [4]: n_col
alpar@9	349
alpar@9	350 -5 number of nonzeros in matrix is negative
alpar@9	351
alpar@9	352 stats [4]: number of nonzeros, p [n_col]
alpar@9	353
alpar@9	354 -6 p [0] is nonzero
alpar@9	355
alpar@9	356 stats [4]: p [0]
alpar@9	357
alpar@9	358 -7 A is too small
alpar@9	359
alpar@9	360 stats [4]: required size
alpar@9	361 stats [5]: actual size (Alen)
alpar@9	362
alpar@9	363 -8 a column has a negative number of entries
alpar@9	364
alpar@9	365 stats [4]: column with < 0 entries
alpar@9	366 stats [5]: number of entries in col
alpar@9	367
alpar@9	368 -9 a row index is out of bounds
alpar@9	369
alpar@9	370 stats [4]: column with bad row index
alpar@9	371 stats [5]: bad row index
alpar@9	372 stats [6]: n_row, # of rows of matrx
alpar@9	373
alpar@9	374 -10 (unused; see symamd.c)
alpar@9	375
alpar@9	376 -999 (unused; see symamd.c)
alpar@9	377
alpar@9	378 Future versions may return more statistics in the stats array.
alpar@9	379
alpar@9	380 Example:
alpar@9	381
alpar@9	382 See http://www.cise.ufl.edu/research/sparse/colamd/example.c
alpar@9	383 for a complete example.
alpar@9	384
alpar@9	385 To order the columns of a 5-by-4 matrix with 11 nonzero entries in
alpar@9	386 the following nonzero pattern
alpar@9	387
alpar@9	388 x 0 x 0
alpar@9	389 x 0 x x
alpar@9	390 0 x x 0
alpar@9	391 0 0 x x
alpar@9	392 x x 0 0
alpar@9	393
alpar@9	394 with default knobs and no output statistics, do the following:
alpar@9	395
alpar@9	396 #include "colamd.h"
alpar@9	397 #define ALEN 100
alpar@9	398 int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
alpar@9	399 int p [ ] = {0, 3, 5, 9, 11} ;
alpar@9	400 int stats [COLAMD_STATS] ;
alpar@9	401 colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
alpar@9	402
alpar@9	403 The permutation is returned in the array p, and A is destroyed.
alpar@9	404
alpar@9	405 ----------------------------------------------------------------------------
alpar@9	406 symamd:
alpar@9	407 ----------------------------------------------------------------------------
alpar@9	408
alpar@9	409 C syntax:
alpar@9	410
alpar@9	411 #include "colamd.h"
alpar@9	412 int symamd (int n, int A, int p, int *perm,
alpar@9	413 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
alpar@9	414 void (allocate) (size_t, size_t), void (release) (void *)) ;
alpar@9	415 UF_long symamd_l (UF_long n, UF_long A, UF_long p, UF_long *perm,
alpar@9	416 double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
alpar@9	417 void (allocate) (size_t, size_t), void (release) (void *)) ;
alpar@9	418
alpar@9	419 Purpose:
alpar@9	420
alpar@9	421 The symamd routine computes an ordering P of a symmetric sparse
alpar@9	422 matrix A such that the Cholesky factorization PAP' = LL' remains
alpar@9	423 sparse. It is based on a column ordering of a matrix M constructed
alpar@9	424 so that the nonzero pattern of M'M is the same as A. The matrix A
alpar@9	425 is assumed to be symmetric; only the strictly lower triangular part
alpar@9	426 is accessed. You must pass your selected memory allocator (usually
alpar@9	427 calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
alpar@9	428 memory for the temporary matrix M.
alpar@9	429
alpar@9	430 Returns:
alpar@9	431
alpar@9	432 TRUE (1) if successful, FALSE (0) otherwise.
alpar@9	433
alpar@9	434 Arguments:
alpar@9	435
alpar@9	436 int n ; Input argument.
alpar@9	437
alpar@9	438 Number of rows and columns in the symmetrix matrix A.
alpar@9	439 Restriction: n >= 0.
alpar@9	440 Symamd returns FALSE if n is negative.
alpar@9	441
alpar@9	442 int A [nnz] ; Input argument.
alpar@9	443
alpar@9	444 A is an integer array of size nnz, where nnz = p [n].
alpar@9	445
alpar@9	446 The row indices of the entries in column c of the matrix are
alpar@9	447 held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a
alpar@9	448 given column c need not be in ascending order, and duplicate
alpar@9	449 row indices may be present. However, symamd will run faster
alpar@9	450 if the columns are in sorted order with no duplicate entries.
alpar@9	451
alpar@9	452 The matrix is 0-based. That is, rows are in the range 0 to
alpar@9	453 n-1, and columns are in the range 0 to n-1. Symamd
alpar@9	454 returns FALSE if any row index is out of range.
alpar@9	455
alpar@9	456 The contents of A are not modified.
alpar@9	457
alpar@9	458 int p [n+1] ; Input argument.
alpar@9	459
alpar@9	460 p is an integer array of size n+1. On input, it holds the
alpar@9	461 "pointers" for the column form of the matrix A. Column c of
alpar@9	462 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
alpar@9	463 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
alpar@9	464 for all c in the range 0 to n-1. The value p [n] is
alpar@9	465 thus the total number of entries in the pattern of the matrix A.
alpar@9	466 Symamd returns FALSE if these conditions are not met.
alpar@9	467
alpar@9	468 The contents of p are not modified.
alpar@9	469
alpar@9	470 int perm [n+1] ; Output argument.
alpar@9	471
alpar@9	472 On output, if symamd returns TRUE, the array perm holds the
alpar@9	473 permutation P, where perm [0] is the first index in the new
alpar@9	474 ordering, and perm [n-1] is the last. That is, perm [k] = j
alpar@9	475 means that row and column j of A is the kth column in PAP',
alpar@9	476 where k is in the range 0 to n-1 (perm [0] = j means
alpar@9	477 that row and column j of A are the first row and column in
alpar@9	478 PAP'). The array is used as a workspace during the ordering,
alpar@9	479 which is why it must be of length n+1, not just n.
alpar@9	480
alpar@9	481 double knobs [COLAMD_KNOBS] ; Input argument.
alpar@9	482
alpar@9	483 See colamd_set_defaults for a description.
alpar@9	484
alpar@9	485 int stats [COLAMD_STATS] ; Output argument.
alpar@9	486
alpar@9	487 Statistics on the ordering, and error status.
alpar@9	488 See colamd.h for related definitions.
alpar@9	489 Symamd returns FALSE if stats is not present.
alpar@9	490
alpar@9	491 stats [0]: number of dense or empty row and columns ignored
alpar@9	492 (and ordered last in the output permutation
alpar@9	493 perm). Note that a row/column can become
alpar@9	494 "empty" if it contains only "dense" and/or
alpar@9	495 "empty" columns/rows.
alpar@9	496
alpar@9	497 stats [1]: (same as stats [0])
alpar@9	498
alpar@9	499 stats [2]: number of garbage collections performed.
alpar@9	500
alpar@9	501 stats [3]: status code. < 0 is an error code.
alpar@9	502 > 1 is a warning or notice.
alpar@9	503
alpar@9	504 0 OK. Each column of the input matrix contained
alpar@9	505 row indices in increasing order, with no
alpar@9	506 duplicates.
alpar@9	507
alpar@9	508 1 OK, but columns of input matrix were jumbled
alpar@9	509 (unsorted columns or duplicate entries). Symamd
alpar@9	510 had to do some extra work to sort the matrix
alpar@9	511 first and remove duplicate entries, but it
alpar@9	512 still was able to return a valid permutation
alpar@9	513 (return value of symamd was TRUE).
alpar@9	514
alpar@9	515 stats [4]: highest numbered column that
alpar@9	516 is unsorted or has duplicate
alpar@9	517 entries.
alpar@9	518 stats [5]: last seen duplicate or
alpar@9	519 unsorted row index.
alpar@9	520 stats [6]: number of duplicate or
alpar@9	521 unsorted row indices.
alpar@9	522
alpar@9	523 -1 A is a null pointer
alpar@9	524
alpar@9	525 -2 p is a null pointer
alpar@9	526
alpar@9	527 -3 (unused, see colamd.c)
alpar@9	528
alpar@9	529 -4 n is negative
alpar@9	530
alpar@9	531 stats [4]: n
alpar@9	532
alpar@9	533 -5 number of nonzeros in matrix is negative
alpar@9	534
alpar@9	535 stats [4]: # of nonzeros (p [n]).
alpar@9	536
alpar@9	537 -6 p [0] is nonzero
alpar@9	538
alpar@9	539 stats [4]: p [0]
alpar@9	540
alpar@9	541 -7 (unused)
alpar@9	542
alpar@9	543 -8 a column has a negative number of entries
alpar@9	544
alpar@9	545 stats [4]: column with < 0 entries
alpar@9	546 stats [5]: number of entries in col
alpar@9	547
alpar@9	548 -9 a row index is out of bounds
alpar@9	549
alpar@9	550 stats [4]: column with bad row index
alpar@9	551 stats [5]: bad row index
alpar@9	552 stats [6]: n_row, # of rows of matrx
alpar@9	553
alpar@9	554 -10 out of memory (unable to allocate temporary
alpar@9	555 workspace for M or count arrays using the
alpar@9	556 "allocate" routine passed into symamd).
alpar@9	557
alpar@9	558 Future versions may return more statistics in the stats array.
alpar@9	559
alpar@9	560 void * (*allocate) (size_t, size_t)
alpar@9	561
alpar@9	562 A pointer to a function providing memory allocation. The
alpar@9	563 allocated memory must be returned initialized to zero. For a
alpar@9	564 C application, this argument should normally be a pointer to
alpar@9	565 calloc. For a MATLAB mexFunction, the routine mxCalloc is
alpar@9	566 passed instead.
alpar@9	567
alpar@9	568 void (*release) (size_t, size_t)
alpar@9	569
alpar@9	570 A pointer to a function that frees memory allocated by the
alpar@9	571 memory allocation routine above. For a C application, this
alpar@9	572 argument should normally be a pointer to free. For a MATLAB
alpar@9	573 mexFunction, the routine mxFree is passed instead.
alpar@9	574
alpar@9	575
alpar@9	576 ----------------------------------------------------------------------------
alpar@9	577 colamd_report:
alpar@9	578 ----------------------------------------------------------------------------
alpar@9	579
alpar@9	580 C syntax:
alpar@9	581
alpar@9	582 #include "colamd.h"
alpar@9	583 colamd_report (int stats [COLAMD_STATS]) ;
alpar@9	584 colamd_l_report (UF_long stats [COLAMD_STATS]) ;
alpar@9	585
alpar@9	586 Purpose:
alpar@9	587
alpar@9	588 Prints the error status and statistics recorded in the stats
alpar@9	589 array on the standard error output (for a standard C routine)
alpar@9	590 or on the MATLAB output (for a mexFunction).
alpar@9	591
alpar@9	592 Arguments:
alpar@9	593
alpar@9	594 int stats [COLAMD_STATS] ; Input only. Statistics from colamd.
alpar@9	595
alpar@9	596
alpar@9	597 ----------------------------------------------------------------------------
alpar@9	598 symamd_report:
alpar@9	599 ----------------------------------------------------------------------------
alpar@9	600
alpar@9	601 C syntax:
alpar@9	602
alpar@9	603 #include "colamd.h"
alpar@9	604 symamd_report (int stats [COLAMD_STATS]) ;
alpar@9	605 symamd_l_report (UF_long stats [COLAMD_STATS]) ;
alpar@9	606
alpar@9	607 Purpose:
alpar@9	608
alpar@9	609 Prints the error status and statistics recorded in the stats
alpar@9	610 array on the standard error output (for a standard C routine)
alpar@9	611 or on the MATLAB output (for a mexFunction).
alpar@9	612
alpar@9	613 Arguments:
alpar@9	614
alpar@9	615 int stats [COLAMD_STATS] ; Input only. Statistics from symamd.
alpar@9	616
alpar@9	617
alpar@9	618 */
alpar@9	619
alpar@9	620 /* ========================================================================== */
alpar@9	621 /* === Scaffolding code definitions ======================================== */
alpar@9	622 /* ========================================================================== */
alpar@9	623
alpar@9	624 /* Ensure that debugging is turned off: */
alpar@9	625 #ifndef NDEBUG
alpar@9	626 #define NDEBUG
alpar@9	627 #endif
alpar@9	628
alpar@9	629 /* turn on debugging by uncommenting the following line
alpar@9	630 #undef NDEBUG
alpar@9	631 */
alpar@9	632
alpar@9	633 /*
alpar@9	634 Our "scaffolding code" philosophy: In our opinion, well-written library
alpar@9	635 code should keep its "debugging" code, and just normally have it turned off
alpar@9	636 by the compiler so as not to interfere with performance. This serves
alpar@9	637 several purposes:
alpar@9	638
alpar@9	639 (1) assertions act as comments to the reader, telling you what the code
alpar@9	640 expects at that point. All assertions will always be true (unless
alpar@9	641 there really is a bug, of course).
alpar@9	642
alpar@9	643 (2) leaving in the scaffolding code assists anyone who would like to modify
alpar@9	644 the code, or understand the algorithm (by reading the debugging output,
alpar@9	645 one can get a glimpse into what the code is doing).
alpar@9	646
alpar@9	647 (3) (gasp!) for actually finding bugs. This code has been heavily tested
alpar@9	648 and "should" be fully functional and bug-free ... but you never know...
alpar@9	649
alpar@9	650 The code will become outrageously slow when debugging is
alpar@9	651 enabled. To control the level of debugging output, set an environment
alpar@9	652 variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging,
alpar@9	653 you should see the following message on the standard output:
alpar@9	654
alpar@9	655 colamd: debug version, D = 1 (THIS WILL BE SLOW!)
alpar@9	656
alpar@9	657 or a similar message for symamd. If you don't, then debugging has not
alpar@9	658 been enabled.
alpar@9	659
alpar@9	660 */
alpar@9	661
alpar@9	662 /* ========================================================================== */
alpar@9	663 /* === Include files ======================================================== */
alpar@9	664 /* ========================================================================== */
alpar@9	665
alpar@9	666 #include "colamd.h"
alpar@9	667
alpar@9	668 #if 0 /* by mao */
alpar@9	669 #include <limits.h>
alpar@9	670 #include <math.h>
alpar@9	671
alpar@9	672 #ifdef MATLAB_MEX_FILE
alpar@9	673 #include "mex.h"
alpar@9	674 #include "matrix.h"
alpar@9	675 #endif /* MATLAB_MEX_FILE */
alpar@9	676
alpar@9	677 #if !defined (NPRINT) \|\| !defined (NDEBUG)
alpar@9	678 #include <stdio.h>
alpar@9	679 #endif
alpar@9	680
alpar@9	681 #ifndef NULL
alpar@9	682 #define NULL ((void *) 0)
alpar@9	683 #endif
alpar@9	684 #endif
alpar@9	685
alpar@9	686 /* ========================================================================== */
alpar@9	687 /* === int or UF_long ======================================================= */
alpar@9	688 /* ========================================================================== */
alpar@9	689
alpar@9	690 #if 0 /* by mao */
alpar@9	691 /* define UF_long */
alpar@9	692 #include "UFconfig.h"
alpar@9	693 #endif
alpar@9	694
alpar@9	695 #ifdef DLONG
alpar@9	696
alpar@9	697 #define Int UF_long
alpar@9	698 #define ID UF_long_id
alpar@9	699 #define Int_MAX UF_long_max
alpar@9	700
alpar@9	701 #define COLAMD_recommended colamd_l_recommended
alpar@9	702 #define COLAMD_set_defaults colamd_l_set_defaults
alpar@9	703 #define COLAMD_MAIN colamd_l
alpar@9	704 #define SYMAMD_MAIN symamd_l
alpar@9	705 #define COLAMD_report colamd_l_report
alpar@9	706 #define SYMAMD_report symamd_l_report
alpar@9	707
alpar@9	708 #else
alpar@9	709
alpar@9	710 #define Int int
alpar@9	711 #define ID "%d"
alpar@9	712 #define Int_MAX INT_MAX
alpar@9	713
alpar@9	714 #define COLAMD_recommended colamd_recommended
alpar@9	715 #define COLAMD_set_defaults colamd_set_defaults
alpar@9	716 #define COLAMD_MAIN colamd
alpar@9	717 #define SYMAMD_MAIN symamd
alpar@9	718 #define COLAMD_report colamd_report
alpar@9	719 #define SYMAMD_report symamd_report
alpar@9	720
alpar@9	721 #endif
alpar@9	722
alpar@9	723 /* ========================================================================== */
alpar@9	724 /* === Row and Column structures ============================================ */
alpar@9	725 /* ========================================================================== */
alpar@9	726
alpar@9	727 /* User code that makes use of the colamd/symamd routines need not directly */
alpar@9	728 /* reference these structures. They are used only for colamd_recommended. */
alpar@9	729
alpar@9	730 typedef struct Colamd_Col_struct
alpar@9	731 {
alpar@9	732 Int start ; /* index for A of first row in this column, or DEAD */
alpar@9	733 /* if column is dead */
alpar@9	734 Int length ; /* number of rows in this column */
alpar@9	735 union
alpar@9	736 {
alpar@9	737 Int thickness ; /* number of original columns represented by this */
alpar@9	738 /* col, if the column is alive */
alpar@9	739 Int parent ; /* parent in parent tree super-column structure, if */
alpar@9	740 /* the column is dead */
alpar@9	741 } shared1 ;
alpar@9	742 union
alpar@9	743 {
alpar@9	744 Int score ; /* the score used to maintain heap, if col is alive */
alpar@9	745 Int order ; /* pivot ordering of this column, if col is dead */
alpar@9	746 } shared2 ;
alpar@9	747 union
alpar@9	748 {
alpar@9	749 Int headhash ; /* head of a hash bucket, if col is at the head of */
alpar@9	750 /* a degree list */
alpar@9	751 Int hash ; /* hash value, if col is not in a degree list */
alpar@9	752 Int prev ; /* previous column in degree list, if col is in a */
alpar@9	753 /* degree list (but not at the head of a degree list) */
alpar@9	754 } shared3 ;
alpar@9	755 union
alpar@9	756 {
alpar@9	757 Int degree_next ; /* next column, if col is in a degree list */
alpar@9	758 Int hash_next ; /* next column, if col is in a hash list */
alpar@9	759 } shared4 ;
alpar@9	760
alpar@9	761 } Colamd_Col ;
alpar@9	762
alpar@9	763 typedef struct Colamd_Row_struct
alpar@9	764 {
alpar@9	765 Int start ; /* index for A of first col in this row */
alpar@9	766 Int length ; /* number of principal columns in this row */
alpar@9	767 union
alpar@9	768 {
alpar@9	769 Int degree ; /* number of principal & non-principal columns in row */
alpar@9	770 Int p ; /* used as a row pointer in init_rows_cols () */
alpar@9	771 } shared1 ;
alpar@9	772 union
alpar@9	773 {
alpar@9	774 Int mark ; /* for computing set differences and marking dead rows*/
alpar@9	775 Int first_column ;/* first column in row (used in garbage collection) */
alpar@9	776 } shared2 ;
alpar@9	777
alpar@9	778 } Colamd_Row ;
alpar@9	779
alpar@9	780 /* ========================================================================== */
alpar@9	781 /* === Definitions ========================================================== */
alpar@9	782 /* ========================================================================== */
alpar@9	783
alpar@9	784 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
alpar@9	785 #define PUBLIC
alpar@9	786 #define PRIVATE static
alpar@9	787
alpar@9	788 #define DENSE_DEGREE(alpha,n) \
alpar@9	789 ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
alpar@9	790
alpar@9	791 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
alpar@9	792 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
alpar@9	793
alpar@9	794 #define ONES_COMPLEMENT(r) (-(r)-1)
alpar@9	795
alpar@9	796 /* -------------------------------------------------------------------------- */
alpar@9	797 /* Change for version 2.1: define TRUE and FALSE only if not yet defined */
alpar@9	798 /* -------------------------------------------------------------------------- */
alpar@9	799
alpar@9	800 #ifndef TRUE
alpar@9	801 #define TRUE (1)
alpar@9	802 #endif
alpar@9	803
alpar@9	804 #ifndef FALSE
alpar@9	805 #define FALSE (0)
alpar@9	806 #endif
alpar@9	807
alpar@9	808 /* -------------------------------------------------------------------------- */
alpar@9	809
alpar@9	810 #define EMPTY (-1)
alpar@9	811
alpar@9	812 /* Row and column status */
alpar@9	813 #define ALIVE (0)
alpar@9	814 #define DEAD (-1)
alpar@9	815
alpar@9	816 /* Column status */
alpar@9	817 #define DEAD_PRINCIPAL (-1)
alpar@9	818 #define DEAD_NON_PRINCIPAL (-2)
alpar@9	819
alpar@9	820 /* Macros for row and column status update and checking. */
alpar@9	821 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
alpar@9	822 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
alpar@9	823 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
alpar@9	824 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
alpar@9	825 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
alpar@9	826 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
alpar@9	827 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
alpar@9	828 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
alpar@9	829 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
alpar@9	830
alpar@9	831 /* ========================================================================== */
alpar@9	832 /* === Colamd reporting mechanism =========================================== */
alpar@9	833 /* ========================================================================== */
alpar@9	834
alpar@9	835 #if defined (MATLAB_MEX_FILE) \|\| defined (MATHWORKS)
alpar@9	836 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
alpar@9	837 #define INDEX(i) ((i)+1)
alpar@9	838 #else
alpar@9	839 /* In C, matrices are 0-based and indices are reported as such in _report /
alpar@9	840 #define INDEX(i) (i)
alpar@9	841 #endif
alpar@9	842
alpar@9	843 /* All output goes through the PRINTF macro. */
alpar@9	844 #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
alpar@9	845
alpar@9	846 /* ========================================================================== */
alpar@9	847 /* === Prototypes of PRIVATE routines ======================================= */
alpar@9	848 /* ========================================================================== */
alpar@9	849
alpar@9	850 PRIVATE Int init_rows_cols
alpar@9	851 (
alpar@9	852 Int n_row,
alpar@9	853 Int n_col,
alpar@9	854 Colamd_Row Row [],
alpar@9	855 Colamd_Col Col [],
alpar@9	856 Int A [],
alpar@9	857 Int p [],
alpar@9	858 Int stats [COLAMD_STATS]
alpar@9	859 ) ;
alpar@9	860
alpar@9	861 PRIVATE void init_scoring
alpar@9	862 (
alpar@9	863 Int n_row,
alpar@9	864 Int n_col,
alpar@9	865 Colamd_Row Row [],
alpar@9	866 Colamd_Col Col [],
alpar@9	867 Int A [],
alpar@9	868 Int head [],
alpar@9	869 double knobs [COLAMD_KNOBS],
alpar@9	870 Int *p_n_row2,
alpar@9	871 Int *p_n_col2,
alpar@9	872 Int *p_max_deg
alpar@9	873 ) ;
alpar@9	874
alpar@9	875 PRIVATE Int find_ordering
alpar@9	876 (
alpar@9	877 Int n_row,
alpar@9	878 Int n_col,
alpar@9	879 Int Alen,
alpar@9	880 Colamd_Row Row [],
alpar@9	881 Colamd_Col Col [],
alpar@9	882 Int A [],
alpar@9	883 Int head [],
alpar@9	884 Int n_col2,
alpar@9	885 Int max_deg,
alpar@9	886 Int pfree,
alpar@9	887 Int aggressive
alpar@9	888 ) ;
alpar@9	889
alpar@9	890 PRIVATE void order_children
alpar@9	891 (
alpar@9	892 Int n_col,
alpar@9	893 Colamd_Col Col [],
alpar@9	894 Int p []
alpar@9	895 ) ;
alpar@9	896
alpar@9	897 PRIVATE void detect_super_cols
alpar@9	898 (
alpar@9	899
alpar@9	900 #ifndef NDEBUG
alpar@9	901 Int n_col,
alpar@9	902 Colamd_Row Row [],
alpar@9	903 #endif /* NDEBUG */
alpar@9	904
alpar@9	905 Colamd_Col Col [],
alpar@9	906 Int A [],
alpar@9	907 Int head [],
alpar@9	908 Int row_start,
alpar@9	909 Int row_length
alpar@9	910 ) ;
alpar@9	911
alpar@9	912 PRIVATE Int garbage_collection
alpar@9	913 (
alpar@9	914 Int n_row,
alpar@9	915 Int n_col,
alpar@9	916 Colamd_Row Row [],
alpar@9	917 Colamd_Col Col [],
alpar@9	918 Int A [],
alpar@9	919 Int *pfree
alpar@9	920 ) ;
alpar@9	921
alpar@9	922 PRIVATE Int clear_mark
alpar@9	923 (
alpar@9	924 Int tag_mark,
alpar@9	925 Int max_mark,
alpar@9	926 Int n_row,
alpar@9	927 Colamd_Row Row []
alpar@9	928 ) ;
alpar@9	929
alpar@9	930 PRIVATE void print_report
alpar@9	931 (
alpar@9	932 char *method,
alpar@9	933 Int stats [COLAMD_STATS]
alpar@9	934 ) ;
alpar@9	935
alpar@9	936 /* ========================================================================== */
alpar@9	937 /* === Debugging prototypes and definitions ================================= */
alpar@9	938 /* ========================================================================== */
alpar@9	939
alpar@9	940 #ifndef NDEBUG
alpar@9	941
alpar@9	942 #if 0 /* by mao */
alpar@9	943 #include <assert.h>
alpar@9	944 #endif
alpar@9	945
alpar@9	946 /* colamd_debug is the ONLY global variable, and is only */
alpar@9	947 /* present when debugging */
alpar@9	948
alpar@9	949 PRIVATE Int colamd_debug = 0 ; /* debug print level */
alpar@9	950
alpar@9	951 #define DEBUG0(params) { PRINTF (params) ; }
alpar@9	952 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
alpar@9	953 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
alpar@9	954 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
alpar@9	955 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
alpar@9	956
alpar@9	957 #if 0 /* by mao */
alpar@9	958 #ifdef MATLAB_MEX_FILE
alpar@9	959 #define ASSERT(expression) (mxAssert ((expression), ""))
alpar@9	960 #else
alpar@9	961 #define ASSERT(expression) (assert (expression))
alpar@9	962 #endif /* MATLAB_MEX_FILE */
alpar@9	963 #else
alpar@9	964 #define ASSERT xassert
alpar@9	965 #endif
alpar@9	966
alpar@9	967 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */
alpar@9	968 (
alpar@9	969 char *method
alpar@9	970 ) ;
alpar@9	971
alpar@9	972 PRIVATE void debug_deg_lists
alpar@9	973 (
alpar@9	974 Int n_row,
alpar@9	975 Int n_col,
alpar@9	976 Colamd_Row Row [],
alpar@9	977 Colamd_Col Col [],
alpar@9	978 Int head [],
alpar@9	979 Int min_score,
alpar@9	980 Int should,
alpar@9	981 Int max_deg
alpar@9	982 ) ;
alpar@9	983
alpar@9	984 PRIVATE void debug_mark
alpar@9	985 (
alpar@9	986 Int n_row,
alpar@9	987 Colamd_Row Row [],
alpar@9	988 Int tag_mark,
alpar@9	989 Int max_mark
alpar@9	990 ) ;
alpar@9	991
alpar@9	992 PRIVATE void debug_matrix
alpar@9	993 (
alpar@9	994 Int n_row,
alpar@9	995 Int n_col,
alpar@9	996 Colamd_Row Row [],
alpar@9	997 Colamd_Col Col [],
alpar@9	998 Int A []
alpar@9	999 ) ;
alpar@9	1000
alpar@9	1001 PRIVATE void debug_structures
alpar@9	1002 (
alpar@9	1003 Int n_row,
alpar@9	1004 Int n_col,
alpar@9	1005 Colamd_Row Row [],
alpar@9	1006 Colamd_Col Col [],
alpar@9	1007 Int A [],
alpar@9	1008 Int n_col2
alpar@9	1009 ) ;
alpar@9	1010
alpar@9	1011 #else /* NDEBUG */
alpar@9	1012
alpar@9	1013 /* === No debugging ========================================================= */
alpar@9	1014
alpar@9	1015 #define DEBUG0(params) ;
alpar@9	1016 #define DEBUG1(params) ;
alpar@9	1017 #define DEBUG2(params) ;
alpar@9	1018 #define DEBUG3(params) ;
alpar@9	1019 #define DEBUG4(params) ;
alpar@9	1020
alpar@9	1021 #define ASSERT(expression)
alpar@9	1022
alpar@9	1023 #endif /* NDEBUG */
alpar@9	1024
alpar@9	1025 /* ========================================================================== */
alpar@9	1026 /* === USER-CALLABLE ROUTINES: ============================================== */
alpar@9	1027 /* ========================================================================== */
alpar@9	1028
alpar@9	1029 /* ========================================================================== */
alpar@9	1030 /* === colamd_recommended =================================================== */
alpar@9	1031 /* ========================================================================== */
alpar@9	1032
alpar@9	1033 /*
alpar@9	1034 The colamd_recommended routine returns the suggested size for Alen. This
alpar@9	1035 value has been determined to provide good balance between the number of
alpar@9	1036 garbage collections and the memory requirements for colamd. If any
alpar@9	1037 argument is negative, or if integer overflow occurs, a 0 is returned as an
alpar@9	1038 error condition. 2*nnz space is required for the row and column
alpar@9	1039 indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
alpar@9	1040 required for the Col and Row arrays, respectively, which are internal to
alpar@9	1041 colamd (roughly 6n_col + 4n_row). An additional n_col space is the
alpar@9	1042 minimal amount of "elbow room", and nnz/5 more space is recommended for
alpar@9	1043 run time efficiency.
alpar@9	1044
alpar@9	1045 Alen is approximately 2.2nnz + 7n_col + 4*n_row + 10.
alpar@9	1046
alpar@9	1047 This function is not needed when using symamd.
alpar@9	1048 */
alpar@9	1049
alpar@9	1050 /* add two values of type size_t, and check for integer overflow */
alpar@9	1051 static size_t t_add (size_t a, size_t b, int *ok)
alpar@9	1052 {
alpar@9	1053 (ok) = (ok) && ((a + b) >= MAX (a,b)) ;
alpar@9	1054 return ((*ok) ? (a + b) : 0) ;
alpar@9	1055 }
alpar@9	1056
alpar@9	1057 /* compute ak where k is a small integer, and check for integer overflow /
alpar@9	1058 static size_t t_mult (size_t a, size_t k, int *ok)
alpar@9	1059 {
alpar@9	1060 size_t i, s = 0 ;
alpar@9	1061 for (i = 0 ; i < k ; i++)
alpar@9	1062 {
alpar@9	1063 s = t_add (s, a, ok) ;
alpar@9	1064 }
alpar@9	1065 return (s) ;
alpar@9	1066 }
alpar@9	1067
alpar@9	1068 /* size of the Col and Row structures */
alpar@9	1069 #define COLAMD_C(n_col,ok) \
alpar@9	1070 ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
alpar@9	1071
alpar@9	1072 #define COLAMD_R(n_row,ok) \
alpar@9	1073 ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
alpar@9	1074
alpar@9	1075
alpar@9	1076 PUBLIC size_t COLAMD_recommended /* returns recommended value of Alen. */
alpar@9	1077 (
alpar@9	1078 /* === Parameters ======================================================= */
alpar@9	1079
alpar@9	1080 Int nnz, /* number of nonzeros in A */
alpar@9	1081 Int n_row, /* number of rows in A */
alpar@9	1082 Int n_col /* number of columns in A */
alpar@9	1083 )
alpar@9	1084 {
alpar@9	1085 size_t s, c, r ;
alpar@9	1086 int ok = TRUE ;
alpar@9	1087 if (nnz < 0 \|\| n_row < 0 \|\| n_col < 0)
alpar@9	1088 {
alpar@9	1089 return (0) ;
alpar@9	1090 }
alpar@9	1091 s = t_mult (nnz, 2, &ok) ; /* 2nnz /
alpar@9	1092 c = COLAMD_C (n_col, &ok) ; /* size of column structures */
alpar@9	1093 r = COLAMD_R (n_row, &ok) ; /* size of row structures */
alpar@9	1094 s = t_add (s, c, &ok) ;
alpar@9	1095 s = t_add (s, r, &ok) ;
alpar@9	1096 s = t_add (s, n_col, &ok) ; /* elbow room */
alpar@9	1097 s = t_add (s, nnz/5, &ok) ; /* elbow room */
alpar@9	1098 ok = ok && (s < Int_MAX) ;
alpar@9	1099 return (ok ? s : 0) ;
alpar@9	1100 }
alpar@9	1101
alpar@9	1102
alpar@9	1103 /* ========================================================================== */
alpar@9	1104 /* === colamd_set_defaults ================================================== */
alpar@9	1105 /* ========================================================================== */
alpar@9	1106
alpar@9	1107 /*
alpar@9	1108 The colamd_set_defaults routine sets the default values of the user-
alpar@9	1109 controllable parameters for colamd and symamd:
alpar@9	1110
alpar@9	1111 Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
alpar@9	1112 entries are removed prior to ordering. Columns with more than
alpar@9	1113 max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
alpar@9	1114 prior to ordering, and placed last in the output column ordering.
alpar@9	1115
alpar@9	1116 Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
alpar@9	1117 entries are removed prior to ordering, and placed last in the
alpar@9	1118 output ordering.
alpar@9	1119
alpar@9	1120 knobs [0] dense row control
alpar@9	1121
alpar@9	1122 knobs [1] dense column control
alpar@9	1123
alpar@9	1124 knobs [2] if nonzero, do aggresive absorption
alpar@9	1125
alpar@9	1126 knobs [3..19] unused, but future versions might use this
alpar@9	1127
alpar@9	1128 */
alpar@9	1129
alpar@9	1130 PUBLIC void COLAMD_set_defaults
alpar@9	1131 (
alpar@9	1132 /* === Parameters ======================================================= */
alpar@9	1133
alpar@9	1134 double knobs [COLAMD_KNOBS] /* knob array */
alpar@9	1135 )
alpar@9	1136 {
alpar@9	1137 /* === Local variables ================================================== */
alpar@9	1138
alpar@9	1139 Int i ;
alpar@9	1140
alpar@9	1141 if (!knobs)
alpar@9	1142 {
alpar@9	1143 return ; /* no knobs to initialize */
alpar@9	1144 }
alpar@9	1145 for (i = 0 ; i < COLAMD_KNOBS ; i++)
alpar@9	1146 {
alpar@9	1147 knobs [i] = 0 ;
alpar@9	1148 }
alpar@9	1149 knobs [COLAMD_DENSE_ROW] = 10 ;
alpar@9	1150 knobs [COLAMD_DENSE_COL] = 10 ;
alpar@9	1151 knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
alpar@9	1152 }
alpar@9	1153
alpar@9	1154
alpar@9	1155 /* ========================================================================== */
alpar@9	1156 /* === symamd =============================================================== */
alpar@9	1157 /* ========================================================================== */
alpar@9	1158
alpar@9	1159 PUBLIC Int SYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */
alpar@9	1160 (
alpar@9	1161 /* === Parameters ======================================================= */
alpar@9	1162
alpar@9	1163 Int n, /* number of rows and columns of A */
alpar@9	1164 Int A [], /* row indices of A */
alpar@9	1165 Int p [], /* column pointers of A */
alpar@9	1166 Int perm [], /* output permutation, size n+1 */
alpar@9	1167 double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */
alpar@9	1168 Int stats [COLAMD_STATS], /* output statistics and error codes */
alpar@9	1169 void * (*allocate) (size_t, size_t),
alpar@9	1170 /* pointer to calloc (ANSI C) or */
alpar@9	1171 /* mxCalloc (for MATLAB mexFunction) */
alpar@9	1172 void (release) (void )
alpar@9	1173 /* pointer to free (ANSI C) or */
alpar@9	1174 /* mxFree (for MATLAB mexFunction) */
alpar@9	1175 )
alpar@9	1176 {
alpar@9	1177 /* === Local variables ================================================== */
alpar@9	1178
alpar@9	1179 Int count ; / length of each column of M, and col pointer*/
alpar@9	1180 Int mark ; / mark array for finding duplicate entries */
alpar@9	1181 Int M ; / row indices of matrix M */
alpar@9	1182 size_t Mlen ; /* length of M */
alpar@9	1183 Int n_row ; /* number of rows in M */
alpar@9	1184 Int nnz ; /* number of entries in A */
alpar@9	1185 Int i ; /* row index of A */
alpar@9	1186 Int j ; /* column index of A */
alpar@9	1187 Int k ; /* row index of M */
alpar@9	1188 Int mnz ; /* number of nonzeros in M */
alpar@9	1189 Int pp ; /* index into a column of A */
alpar@9	1190 Int last_row ; /* last row seen in the current column */
alpar@9	1191 Int length ; /* number of nonzeros in a column */
alpar@9	1192
alpar@9	1193 double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */
alpar@9	1194 double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */
alpar@9	1195
alpar@9	1196 #ifndef NDEBUG
alpar@9	1197 colamd_get_debug ("symamd") ;
alpar@9	1198 #endif /* NDEBUG */
alpar@9	1199
alpar@9	1200 /* === Check the input arguments ======================================== */
alpar@9	1201
alpar@9	1202 if (!stats)
alpar@9	1203 {
alpar@9	1204 DEBUG0 (("symamd: stats not present\n")) ;
alpar@9	1205 return (FALSE) ;
alpar@9	1206 }
alpar@9	1207 for (i = 0 ; i < COLAMD_STATS ; i++)
alpar@9	1208 {
alpar@9	1209 stats [i] = 0 ;
alpar@9	1210 }
alpar@9	1211 stats [COLAMD_STATUS] = COLAMD_OK ;
alpar@9	1212 stats [COLAMD_INFO1] = -1 ;
alpar@9	1213 stats [COLAMD_INFO2] = -1 ;
alpar@9	1214
alpar@9	1215 if (!A)
alpar@9	1216 {
alpar@9	1217 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
alpar@9	1218 DEBUG0 (("symamd: A not present\n")) ;
alpar@9	1219 return (FALSE) ;
alpar@9	1220 }
alpar@9	1221
alpar@9	1222 if (!p) /* p is not present */
alpar@9	1223 {
alpar@9	1224 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
alpar@9	1225 DEBUG0 (("symamd: p not present\n")) ;
alpar@9	1226 return (FALSE) ;
alpar@9	1227 }
alpar@9	1228
alpar@9	1229 if (n < 0) /* n must be >= 0 */
alpar@9	1230 {
alpar@9	1231 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
alpar@9	1232 stats [COLAMD_INFO1] = n ;
alpar@9	1233 DEBUG0 (("symamd: n negative %d\n", n)) ;
alpar@9	1234 return (FALSE) ;
alpar@9	1235 }
alpar@9	1236
alpar@9	1237 nnz = p [n] ;
alpar@9	1238 if (nnz < 0) /* nnz must be >= 0 */
alpar@9	1239 {
alpar@9	1240 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
alpar@9	1241 stats [COLAMD_INFO1] = nnz ;
alpar@9	1242 DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
alpar@9	1243 return (FALSE) ;
alpar@9	1244 }
alpar@9	1245
alpar@9	1246 if (p [0] != 0)
alpar@9	1247 {
alpar@9	1248 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
alpar@9	1249 stats [COLAMD_INFO1] = p [0] ;
alpar@9	1250 DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
alpar@9	1251 return (FALSE) ;
alpar@9	1252 }
alpar@9	1253
alpar@9	1254 /* === If no knobs, set default knobs =================================== */
alpar@9	1255
alpar@9	1256 if (!knobs)
alpar@9	1257 {
alpar@9	1258 COLAMD_set_defaults (default_knobs) ;
alpar@9	1259 knobs = default_knobs ;
alpar@9	1260 }
alpar@9	1261
alpar@9	1262 /* === Allocate count and mark ========================================== */
alpar@9	1263
alpar@9	1264 count = (Int ) ((allocate) (n+1, sizeof (Int))) ;
alpar@9	1265 if (!count)
alpar@9	1266 {
alpar@9	1267 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
alpar@9	1268 DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
alpar@9	1269 return (FALSE) ;
alpar@9	1270 }
alpar@9	1271
alpar@9	1272 mark = (Int ) ((allocate) (n+1, sizeof (Int))) ;
alpar@9	1273 if (!mark)
alpar@9	1274 {
alpar@9	1275 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
alpar@9	1276 (release) ((void ) count) ;
alpar@9	1277 DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
alpar@9	1278 return (FALSE) ;
alpar@9	1279 }
alpar@9	1280
alpar@9	1281 /* === Compute column counts of M, check if A is valid ================== */
alpar@9	1282
alpar@9	1283 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
alpar@9	1284
alpar@9	1285 for (i = 0 ; i < n ; i++)
alpar@9	1286 {
alpar@9	1287 mark [i] = -1 ;
alpar@9	1288 }
alpar@9	1289
alpar@9	1290 for (j = 0 ; j < n ; j++)
alpar@9	1291 {
alpar@9	1292 last_row = -1 ;
alpar@9	1293
alpar@9	1294 length = p [j+1] - p [j] ;
alpar@9	1295 if (length < 0)
alpar@9	1296 {
alpar@9	1297 /* column pointers must be non-decreasing */
alpar@9	1298 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
alpar@9	1299 stats [COLAMD_INFO1] = j ;
alpar@9	1300 stats [COLAMD_INFO2] = length ;
alpar@9	1301 (release) ((void ) count) ;
alpar@9	1302 (release) ((void ) mark) ;
alpar@9	1303 DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
alpar@9	1304 return (FALSE) ;
alpar@9	1305 }
alpar@9	1306
alpar@9	1307 for (pp = p [j] ; pp < p [j+1] ; pp++)
alpar@9	1308 {
alpar@9	1309 i = A [pp] ;
alpar@9	1310 if (i < 0 \|\| i >= n)
alpar@9	1311 {
alpar@9	1312 /* row index i, in column j, is out of bounds */
alpar@9	1313 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
alpar@9	1314 stats [COLAMD_INFO1] = j ;
alpar@9	1315 stats [COLAMD_INFO2] = i ;
alpar@9	1316 stats [COLAMD_INFO3] = n ;
alpar@9	1317 (release) ((void ) count) ;
alpar@9	1318 (release) ((void ) mark) ;
alpar@9	1319 DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
alpar@9	1320 return (FALSE) ;
alpar@9	1321 }
alpar@9	1322
alpar@9	1323 if (i <= last_row \|\| mark [i] == j)
alpar@9	1324 {
alpar@9	1325 /* row index is unsorted or repeated (or both), thus col */
alpar@9	1326 /* is jumbled. This is a notice, not an error condition. */
alpar@9	1327 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
alpar@9	1328 stats [COLAMD_INFO1] = j ;
alpar@9	1329 stats [COLAMD_INFO2] = i ;
alpar@9	1330 (stats [COLAMD_INFO3]) ++ ;
alpar@9	1331 DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
alpar@9	1332 }
alpar@9	1333
alpar@9	1334 if (i > j && mark [i] != j)
alpar@9	1335 {
alpar@9	1336 /* row k of M will contain column indices i and j */
alpar@9	1337 count [i]++ ;
alpar@9	1338 count [j]++ ;
alpar@9	1339 }
alpar@9	1340
alpar@9	1341 /* mark the row as having been seen in this column */
alpar@9	1342 mark [i] = j ;
alpar@9	1343
alpar@9	1344 last_row = i ;
alpar@9	1345 }
alpar@9	1346 }
alpar@9	1347
alpar@9	1348 /* v2.4: removed free(mark) */
alpar@9	1349
alpar@9	1350 /* === Compute column pointers of M ===================================== */
alpar@9	1351
alpar@9	1352 /* use output permutation, perm, for column pointers of M */
alpar@9	1353 perm [0] = 0 ;
alpar@9	1354 for (j = 1 ; j <= n ; j++)
alpar@9	1355 {
alpar@9	1356 perm [j] = perm [j-1] + count [j-1] ;
alpar@9	1357 }
alpar@9	1358 for (j = 0 ; j < n ; j++)
alpar@9	1359 {
alpar@9	1360 count [j] = perm [j] ;
alpar@9	1361 }
alpar@9	1362
alpar@9	1363 /* === Construct M ====================================================== */
alpar@9	1364
alpar@9	1365 mnz = perm [n] ;
alpar@9	1366 n_row = mnz / 2 ;
alpar@9	1367 Mlen = COLAMD_recommended (mnz, n_row, n) ;
alpar@9	1368 M = (Int ) ((allocate) (Mlen, sizeof (Int))) ;
alpar@9	1369 DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
alpar@9	1370 n_row, n, mnz, (double) Mlen)) ;
alpar@9	1371
alpar@9	1372 if (!M)
alpar@9	1373 {
alpar@9	1374 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
alpar@9	1375 (release) ((void ) count) ;
alpar@9	1376 (release) ((void ) mark) ;
alpar@9	1377 DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
alpar@9	1378 return (FALSE) ;
alpar@9	1379 }
alpar@9	1380
alpar@9	1381 k = 0 ;
alpar@9	1382
alpar@9	1383 if (stats [COLAMD_STATUS] == COLAMD_OK)
alpar@9	1384 {
alpar@9	1385 /* Matrix is OK */
alpar@9	1386 for (j = 0 ; j < n ; j++)
alpar@9	1387 {
alpar@9	1388 ASSERT (p [j+1] - p [j] >= 0) ;
alpar@9	1389 for (pp = p [j] ; pp < p [j+1] ; pp++)
alpar@9	1390 {
alpar@9	1391 i = A [pp] ;
alpar@9	1392 ASSERT (i >= 0 && i < n) ;
alpar@9	1393 if (i > j)
alpar@9	1394 {
alpar@9	1395 /* row k of M contains column indices i and j */
alpar@9	1396 M [count [i]++] = k ;
alpar@9	1397 M [count [j]++] = k ;
alpar@9	1398 k++ ;
alpar@9	1399 }
alpar@9	1400 }
alpar@9	1401 }
alpar@9	1402 }
alpar@9	1403 else
alpar@9	1404 {
alpar@9	1405 /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */
alpar@9	1406 DEBUG0 (("symamd: Duplicates in A.\n")) ;
alpar@9	1407 for (i = 0 ; i < n ; i++)
alpar@9	1408 {
alpar@9	1409 mark [i] = -1 ;
alpar@9	1410 }
alpar@9	1411 for (j = 0 ; j < n ; j++)
alpar@9	1412 {
alpar@9	1413 ASSERT (p [j+1] - p [j] >= 0) ;
alpar@9	1414 for (pp = p [j] ; pp < p [j+1] ; pp++)
alpar@9	1415 {
alpar@9	1416 i = A [pp] ;
alpar@9	1417 ASSERT (i >= 0 && i < n) ;
alpar@9	1418 if (i > j && mark [i] != j)
alpar@9	1419 {
alpar@9	1420 /* row k of M contains column indices i and j */
alpar@9	1421 M [count [i]++] = k ;
alpar@9	1422 M [count [j]++] = k ;
alpar@9	1423 k++ ;
alpar@9	1424 mark [i] = j ;
alpar@9	1425 }
alpar@9	1426 }
alpar@9	1427 }
alpar@9	1428 /* v2.4: free(mark) moved below */
alpar@9	1429 }
alpar@9	1430
alpar@9	1431 /* count and mark no longer needed */
alpar@9	1432 (release) ((void ) count) ;
alpar@9	1433 (release) ((void ) mark) ; /* v2.4: free (mark) moved here */
alpar@9	1434 ASSERT (k == n_row) ;
alpar@9	1435
alpar@9	1436 /* === Adjust the knobs for M =========================================== */
alpar@9	1437
alpar@9	1438 for (i = 0 ; i < COLAMD_KNOBS ; i++)
alpar@9	1439 {
alpar@9	1440 cknobs [i] = knobs [i] ;
alpar@9	1441 }
alpar@9	1442
alpar@9	1443 /* there are no dense rows in M */
alpar@9	1444 cknobs [COLAMD_DENSE_ROW] = -1 ;
alpar@9	1445 cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
alpar@9	1446
alpar@9	1447 /* === Order the columns of M =========================================== */
alpar@9	1448
alpar@9	1449 /* v2.4: colamd cannot fail here, so the error check is removed */
alpar@9	1450 (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
alpar@9	1451
alpar@9	1452 /* Note that the output permutation is now in perm */
alpar@9	1453
alpar@9	1454 /* === get the statistics for symamd from colamd ======================== */
alpar@9	1455
alpar@9	1456 /* a dense column in colamd means a dense row and col in symamd */
alpar@9	1457 stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
alpar@9	1458
alpar@9	1459 /* === Free M =========================================================== */
alpar@9	1460
alpar@9	1461 (release) ((void ) M) ;
alpar@9	1462 DEBUG0 (("symamd: done.\n")) ;
alpar@9	1463 return (TRUE) ;
alpar@9	1464
alpar@9	1465 }
alpar@9	1466
alpar@9	1467 /* ========================================================================== */
alpar@9	1468 /* === colamd =============================================================== */
alpar@9	1469 /* ========================================================================== */
alpar@9	1470
alpar@9	1471 /*
alpar@9	1472 The colamd routine computes a column ordering Q of a sparse matrix
alpar@9	1473 A such that the LU factorization P(AQ) = LU remains sparse, where P is
alpar@9	1474 selected via partial pivoting. The routine can also be viewed as
alpar@9	1475 providing a permutation Q such that the Cholesky factorization
alpar@9	1476 (AQ)'(AQ) = LL' remains sparse.
alpar@9	1477 */
alpar@9	1478
alpar@9	1479 PUBLIC Int COLAMD_MAIN /* returns TRUE if successful, FALSE otherwise*/
alpar@9	1480 (
alpar@9	1481 /* === Parameters ======================================================= */
alpar@9	1482
alpar@9	1483 Int n_row, /* number of rows in A */
alpar@9	1484 Int n_col, /* number of columns in A */
alpar@9	1485 Int Alen, /* length of A */
alpar@9	1486 Int A [], /* row indices of A */
alpar@9	1487 Int p [], /* pointers to columns in A */
alpar@9	1488 double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
alpar@9	1489 Int stats [COLAMD_STATS] /* output statistics and error codes */
alpar@9	1490 )
alpar@9	1491 {
alpar@9	1492 /* === Local variables ================================================== */
alpar@9	1493
alpar@9	1494 Int i ; /* loop index */
alpar@9	1495 Int nnz ; /* nonzeros in A */
alpar@9	1496 size_t Row_size ; /* size of Row [], in integers */
alpar@9	1497 size_t Col_size ; /* size of Col [], in integers */
alpar@9	1498 size_t need ; /* minimum required length of A */
alpar@9	1499 Colamd_Row Row ; / pointer into A of Row [0..n_row] array */
alpar@9	1500 Colamd_Col Col ; / pointer into A of Col [0..n_col] array */
alpar@9	1501 Int n_col2 ; /* number of non-dense, non-empty columns */
alpar@9	1502 Int n_row2 ; /* number of non-dense, non-empty rows */
alpar@9	1503 Int ngarbage ; /* number of garbage collections performed */
alpar@9	1504 Int max_deg ; /* maximum row degree */
alpar@9	1505 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
alpar@9	1506 Int aggressive ; /* do aggressive absorption */
alpar@9	1507 int ok ;
alpar@9	1508
alpar@9	1509 #ifndef NDEBUG
alpar@9	1510 colamd_get_debug ("colamd") ;
alpar@9	1511 #endif /* NDEBUG */
alpar@9	1512
alpar@9	1513 /* === Check the input arguments ======================================== */
alpar@9	1514
alpar@9	1515 if (!stats)
alpar@9	1516 {
alpar@9	1517 DEBUG0 (("colamd: stats not present\n")) ;
alpar@9	1518 return (FALSE) ;
alpar@9	1519 }
alpar@9	1520 for (i = 0 ; i < COLAMD_STATS ; i++)
alpar@9	1521 {
alpar@9	1522 stats [i] = 0 ;
alpar@9	1523 }
alpar@9	1524 stats [COLAMD_STATUS] = COLAMD_OK ;
alpar@9	1525 stats [COLAMD_INFO1] = -1 ;
alpar@9	1526 stats [COLAMD_INFO2] = -1 ;
alpar@9	1527
alpar@9	1528 if (!A) /* A is not present */
alpar@9	1529 {
alpar@9	1530 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
alpar@9	1531 DEBUG0 (("colamd: A not present\n")) ;
alpar@9	1532 return (FALSE) ;
alpar@9	1533 }
alpar@9	1534
alpar@9	1535 if (!p) /* p is not present */
alpar@9	1536 {
alpar@9	1537 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
alpar@9	1538 DEBUG0 (("colamd: p not present\n")) ;
alpar@9	1539 return (FALSE) ;
alpar@9	1540 }
alpar@9	1541
alpar@9	1542 if (n_row < 0) /* n_row must be >= 0 */
alpar@9	1543 {
alpar@9	1544 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
alpar@9	1545 stats [COLAMD_INFO1] = n_row ;
alpar@9	1546 DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
alpar@9	1547 return (FALSE) ;
alpar@9	1548 }
alpar@9	1549
alpar@9	1550 if (n_col < 0) /* n_col must be >= 0 */
alpar@9	1551 {
alpar@9	1552 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
alpar@9	1553 stats [COLAMD_INFO1] = n_col ;
alpar@9	1554 DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
alpar@9	1555 return (FALSE) ;
alpar@9	1556 }
alpar@9	1557
alpar@9	1558 nnz = p [n_col] ;
alpar@9	1559 if (nnz < 0) /* nnz must be >= 0 */
alpar@9	1560 {
alpar@9	1561 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
alpar@9	1562 stats [COLAMD_INFO1] = nnz ;
alpar@9	1563 DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
alpar@9	1564 return (FALSE) ;
alpar@9	1565 }
alpar@9	1566
alpar@9	1567 if (p [0] != 0)
alpar@9	1568 {
alpar@9	1569 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
alpar@9	1570 stats [COLAMD_INFO1] = p [0] ;
alpar@9	1571 DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
alpar@9	1572 return (FALSE) ;
alpar@9	1573 }
alpar@9	1574
alpar@9	1575 /* === If no knobs, set default knobs =================================== */
alpar@9	1576
alpar@9	1577 if (!knobs)
alpar@9	1578 {
alpar@9	1579 COLAMD_set_defaults (default_knobs) ;
alpar@9	1580 knobs = default_knobs ;
alpar@9	1581 }
alpar@9	1582
alpar@9	1583 aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
alpar@9	1584
alpar@9	1585 /* === Allocate the Row and Col arrays from array A ===================== */
alpar@9	1586
alpar@9	1587 ok = TRUE ;
alpar@9	1588 Col_size = COLAMD_C (n_col, &ok) ; /* size of Col array of structs */
alpar@9	1589 Row_size = COLAMD_R (n_row, &ok) ; /* size of Row array of structs */
alpar@9	1590
alpar@9	1591 /* need = 2nnz + n_col + Col_size + Row_size ; /
alpar@9	1592 need = t_mult (nnz, 2, &ok) ;
alpar@9	1593 need = t_add (need, n_col, &ok) ;
alpar@9	1594 need = t_add (need, Col_size, &ok) ;
alpar@9	1595 need = t_add (need, Row_size, &ok) ;
alpar@9	1596
alpar@9	1597 if (!ok \|\| need > (size_t) Alen \|\| need > Int_MAX)
alpar@9	1598 {
alpar@9	1599 /* not enough space in array A to perform the ordering */
alpar@9	1600 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
alpar@9	1601 stats [COLAMD_INFO1] = need ;
alpar@9	1602 stats [COLAMD_INFO2] = Alen ;
alpar@9	1603 DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
alpar@9	1604 return (FALSE) ;
alpar@9	1605 }
alpar@9	1606
alpar@9	1607 Alen -= Col_size + Row_size ;
alpar@9	1608 Col = (Colamd_Col *) &A [Alen] ;
alpar@9	1609 Row = (Colamd_Row *) &A [Alen + Col_size] ;
alpar@9	1610
alpar@9	1611 /* === Construct the row and column data structures ===================== */
alpar@9	1612
alpar@9	1613 if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
alpar@9	1614 {
alpar@9	1615 /* input matrix is invalid */
alpar@9	1616 DEBUG0 (("colamd: Matrix invalid\n")) ;
alpar@9	1617 return (FALSE) ;
alpar@9	1618 }
alpar@9	1619
alpar@9	1620 /* === Initialize scores, kill dense rows/columns ======================= */
alpar@9	1621
alpar@9	1622 init_scoring (n_row, n_col, Row, Col, A, p, knobs,
alpar@9	1623 &n_row2, &n_col2, &max_deg) ;
alpar@9	1624
alpar@9	1625 /* === Order the supercolumns =========================================== */
alpar@9	1626
alpar@9	1627 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
alpar@9	1628 n_col2, max_deg, 2*nnz, aggressive) ;
alpar@9	1629
alpar@9	1630 /* === Order the non-principal columns ================================== */
alpar@9	1631
alpar@9	1632 order_children (n_col, Col, p) ;
alpar@9	1633
alpar@9	1634 /* === Return statistics in stats ======================================= */
alpar@9	1635
alpar@9	1636 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
alpar@9	1637 stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
alpar@9	1638 stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
alpar@9	1639 DEBUG0 (("colamd: done.\n")) ;
alpar@9	1640 return (TRUE) ;
alpar@9	1641 }
alpar@9	1642
alpar@9	1643
alpar@9	1644 /* ========================================================================== */
alpar@9	1645 /* === colamd_report ======================================================== */
alpar@9	1646 /* ========================================================================== */
alpar@9	1647
alpar@9	1648 PUBLIC void COLAMD_report
alpar@9	1649 (
alpar@9	1650 Int stats [COLAMD_STATS]
alpar@9	1651 )
alpar@9	1652 {
alpar@9	1653 print_report ("colamd", stats) ;
alpar@9	1654 }
alpar@9	1655
alpar@9	1656
alpar@9	1657 /* ========================================================================== */
alpar@9	1658 /* === symamd_report ======================================================== */
alpar@9	1659 /* ========================================================================== */
alpar@9	1660
alpar@9	1661 PUBLIC void SYMAMD_report
alpar@9	1662 (
alpar@9	1663 Int stats [COLAMD_STATS]
alpar@9	1664 )
alpar@9	1665 {
alpar@9	1666 print_report ("symamd", stats) ;
alpar@9	1667 }
alpar@9	1668
alpar@9	1669
alpar@9	1670
alpar@9	1671 /* ========================================================================== */
alpar@9	1672 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
alpar@9	1673 /* ========================================================================== */
alpar@9	1674
alpar@9	1675 /* There are no user-callable routines beyond this point in the file */
alpar@9	1676
alpar@9	1677
alpar@9	1678 /* ========================================================================== */
alpar@9	1679 /* === init_rows_cols ======================================================= */
alpar@9	1680 /* ========================================================================== */
alpar@9	1681
alpar@9	1682 /*
alpar@9	1683 Takes the column form of the matrix in A and creates the row form of the
alpar@9	1684 matrix. Also, row and column attributes are stored in the Col and Row
alpar@9	1685 structs. If the columns are un-sorted or contain duplicate row indices,
alpar@9	1686 this routine will also sort and remove duplicate row indices from the
alpar@9	1687 column form of the matrix. Returns FALSE if the matrix is invalid,
alpar@9	1688 TRUE otherwise. Not user-callable.
alpar@9	1689 */
alpar@9	1690
alpar@9	1691 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
alpar@9	1692 (
alpar@9	1693 /* === Parameters ======================================================= */
alpar@9	1694
alpar@9	1695 Int n_row, /* number of rows of A */
alpar@9	1696 Int n_col, /* number of columns of A */
alpar@9	1697 Colamd_Row Row [], /* of size n_row+1 */
alpar@9	1698 Colamd_Col Col [], /* of size n_col+1 */
alpar@9	1699 Int A [], /* row indices of A, of size Alen */
alpar@9	1700 Int p [], /* pointers to columns in A, of size n_col+1 */
alpar@9	1701 Int stats [COLAMD_STATS] /* colamd statistics */
alpar@9	1702 )
alpar@9	1703 {
alpar@9	1704 /* === Local variables ================================================== */
alpar@9	1705
alpar@9	1706 Int col ; /* a column index */
alpar@9	1707 Int row ; /* a row index */
alpar@9	1708 Int cp ; / a column pointer */
alpar@9	1709 Int cp_end ; / a pointer to the end of a column */
alpar@9	1710 Int rp ; / a row pointer */
alpar@9	1711 Int rp_end ; / a pointer to the end of a row */
alpar@9	1712 Int last_row ; /* previous row */
alpar@9	1713
alpar@9	1714 /* === Initialize columns, and check column pointers ==================== */
alpar@9	1715
alpar@9	1716 for (col = 0 ; col < n_col ; col++)
alpar@9	1717 {
alpar@9	1718 Col [col].start = p [col] ;
alpar@9	1719 Col [col].length = p [col+1] - p [col] ;
alpar@9	1720
alpar@9	1721 if (Col [col].length < 0)
alpar@9	1722 {
alpar@9	1723 /* column pointers must be non-decreasing */
alpar@9	1724 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
alpar@9	1725 stats [COLAMD_INFO1] = col ;
alpar@9	1726 stats [COLAMD_INFO2] = Col [col].length ;
alpar@9	1727 DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
alpar@9	1728 return (FALSE) ;
alpar@9	1729 }
alpar@9	1730
alpar@9	1731 Col [col].shared1.thickness = 1 ;
alpar@9	1732 Col [col].shared2.score = 0 ;
alpar@9	1733 Col [col].shared3.prev = EMPTY ;
alpar@9	1734 Col [col].shared4.degree_next = EMPTY ;
alpar@9	1735 }
alpar@9	1736
alpar@9	1737 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
alpar@9	1738
alpar@9	1739 /* === Scan columns, compute row degrees, and check row indices ========= */
alpar@9	1740
alpar@9	1741 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
alpar@9	1742
alpar@9	1743 for (row = 0 ; row < n_row ; row++)
alpar@9	1744 {
alpar@9	1745 Row [row].length = 0 ;
alpar@9	1746 Row [row].shared2.mark = -1 ;
alpar@9	1747 }
alpar@9	1748
alpar@9	1749 for (col = 0 ; col < n_col ; col++)
alpar@9	1750 {
alpar@9	1751 last_row = -1 ;
alpar@9	1752
alpar@9	1753 cp = &A [p [col]] ;
alpar@9	1754 cp_end = &A [p [col+1]] ;
alpar@9	1755
alpar@9	1756 while (cp < cp_end)
alpar@9	1757 {
alpar@9	1758 row = *cp++ ;
alpar@9	1759
alpar@9	1760 /* make sure row indices within range */
alpar@9	1761 if (row < 0 \|\| row >= n_row)
alpar@9	1762 {
alpar@9	1763 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
alpar@9	1764 stats [COLAMD_INFO1] = col ;
alpar@9	1765 stats [COLAMD_INFO2] = row ;
alpar@9	1766 stats [COLAMD_INFO3] = n_row ;
alpar@9	1767 DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
alpar@9	1768 return (FALSE) ;
alpar@9	1769 }
alpar@9	1770
alpar@9	1771 if (row <= last_row \|\| Row [row].shared2.mark == col)
alpar@9	1772 {
alpar@9	1773 /* row index are unsorted or repeated (or both), thus col */
alpar@9	1774 /* is jumbled. This is a notice, not an error condition. */
alpar@9	1775 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
alpar@9	1776 stats [COLAMD_INFO1] = col ;
alpar@9	1777 stats [COLAMD_INFO2] = row ;
alpar@9	1778 (stats [COLAMD_INFO3]) ++ ;
alpar@9	1779 DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
alpar@9	1780 }
alpar@9	1781
alpar@9	1782 if (Row [row].shared2.mark != col)
alpar@9	1783 {
alpar@9	1784 Row [row].length++ ;
alpar@9	1785 }
alpar@9	1786 else
alpar@9	1787 {
alpar@9	1788 /* this is a repeated entry in the column, */
alpar@9	1789 /* it will be removed */
alpar@9	1790 Col [col].length-- ;
alpar@9	1791 }
alpar@9	1792
alpar@9	1793 /* mark the row as having been seen in this column */
alpar@9	1794 Row [row].shared2.mark = col ;
alpar@9	1795
alpar@9	1796 last_row = row ;
alpar@9	1797 }
alpar@9	1798 }
alpar@9	1799
alpar@9	1800 /* === Compute row pointers ============================================= */
alpar@9	1801
alpar@9	1802 /* row form of the matrix starts directly after the column */
alpar@9	1803 /* form of matrix in A */
alpar@9	1804 Row [0].start = p [n_col] ;
alpar@9	1805 Row [0].shared1.p = Row [0].start ;
alpar@9	1806 Row [0].shared2.mark = -1 ;
alpar@9	1807 for (row = 1 ; row < n_row ; row++)
alpar@9	1808 {
alpar@9	1809 Row [row].start = Row [row-1].start + Row [row-1].length ;
alpar@9	1810 Row [row].shared1.p = Row [row].start ;
alpar@9	1811 Row [row].shared2.mark = -1 ;
alpar@9	1812 }
alpar@9	1813
alpar@9	1814 /* === Create row form ================================================== */
alpar@9	1815
alpar@9	1816 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
alpar@9	1817 {
alpar@9	1818 /* if cols jumbled, watch for repeated row indices */
alpar@9	1819 for (col = 0 ; col < n_col ; col++)
alpar@9	1820 {
alpar@9	1821 cp = &A [p [col]] ;
alpar@9	1822 cp_end = &A [p [col+1]] ;
alpar@9	1823 while (cp < cp_end)
alpar@9	1824 {
alpar@9	1825 row = *cp++ ;
alpar@9	1826 if (Row [row].shared2.mark != col)
alpar@9	1827 {
alpar@9	1828 A [(Row [row].shared1.p)++] = col ;
alpar@9	1829 Row [row].shared2.mark = col ;
alpar@9	1830 }
alpar@9	1831 }
alpar@9	1832 }
alpar@9	1833 }
alpar@9	1834 else
alpar@9	1835 {
alpar@9	1836 /* if cols not jumbled, we don't need the mark (this is faster) */
alpar@9	1837 for (col = 0 ; col < n_col ; col++)
alpar@9	1838 {
alpar@9	1839 cp = &A [p [col]] ;
alpar@9	1840 cp_end = &A [p [col+1]] ;
alpar@9	1841 while (cp < cp_end)
alpar@9	1842 {
alpar@9	1843 A [(Row [*cp++].shared1.p)++] = col ;
alpar@9	1844 }
alpar@9	1845 }
alpar@9	1846 }
alpar@9	1847
alpar@9	1848 /* === Clear the row marks and set row degrees ========================== */
alpar@9	1849
alpar@9	1850 for (row = 0 ; row < n_row ; row++)
alpar@9	1851 {
alpar@9	1852 Row [row].shared2.mark = 0 ;
alpar@9	1853 Row [row].shared1.degree = Row [row].length ;
alpar@9	1854 }
alpar@9	1855
alpar@9	1856 /* === See if we need to re-create columns ============================== */
alpar@9	1857
alpar@9	1858 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
alpar@9	1859 {
alpar@9	1860 DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
alpar@9	1861
alpar@9	1862 #ifndef NDEBUG
alpar@9	1863 /* make sure column lengths are correct */
alpar@9	1864 for (col = 0 ; col < n_col ; col++)
alpar@9	1865 {
alpar@9	1866 p [col] = Col [col].length ;
alpar@9	1867 }
alpar@9	1868 for (row = 0 ; row < n_row ; row++)
alpar@9	1869 {
alpar@9	1870 rp = &A [Row [row].start] ;
alpar@9	1871 rp_end = rp + Row [row].length ;
alpar@9	1872 while (rp < rp_end)
alpar@9	1873 {
alpar@9	1874 p [*rp++]-- ;
alpar@9	1875 }
alpar@9	1876 }
alpar@9	1877 for (col = 0 ; col < n_col ; col++)
alpar@9	1878 {
alpar@9	1879 ASSERT (p [col] == 0) ;
alpar@9	1880 }
alpar@9	1881 /* now p is all zero (different than when debugging is turned off) */
alpar@9	1882 #endif /* NDEBUG */
alpar@9	1883
alpar@9	1884 /* === Compute col pointers ========================================= */
alpar@9	1885
alpar@9	1886 /* col form of the matrix starts at A [0]. */
alpar@9	1887 /* Note, we may have a gap between the col form and the row */
alpar@9	1888 /* form if there were duplicate entries, if so, it will be */
alpar@9	1889 /* removed upon the first garbage collection */
alpar@9	1890 Col [0].start = 0 ;
alpar@9	1891 p [0] = Col [0].start ;
alpar@9	1892 for (col = 1 ; col < n_col ; col++)
alpar@9	1893 {
alpar@9	1894 /* note that the lengths here are for pruned columns, i.e. */
alpar@9	1895 /* no duplicate row indices will exist for these columns */
alpar@9	1896 Col [col].start = Col [col-1].start + Col [col-1].length ;
alpar@9	1897 p [col] = Col [col].start ;
alpar@9	1898 }
alpar@9	1899
alpar@9	1900 /* === Re-create col form =========================================== */
alpar@9	1901
alpar@9	1902 for (row = 0 ; row < n_row ; row++)
alpar@9	1903 {
alpar@9	1904 rp = &A [Row [row].start] ;
alpar@9	1905 rp_end = rp + Row [row].length ;
alpar@9	1906 while (rp < rp_end)
alpar@9	1907 {
alpar@9	1908 A [(p [*rp++])++] = row ;
alpar@9	1909 }
alpar@9	1910 }
alpar@9	1911 }
alpar@9	1912
alpar@9	1913 /* === Done. Matrix is not (or no longer) jumbled ====================== */
alpar@9	1914
alpar@9	1915 return (TRUE) ;
alpar@9	1916 }
alpar@9	1917
alpar@9	1918
alpar@9	1919 /* ========================================================================== */
alpar@9	1920 /* === init_scoring ========================================================= */
alpar@9	1921 /* ========================================================================== */
alpar@9	1922
alpar@9	1923 /*
alpar@9	1924 Kills dense or empty columns and rows, calculates an initial score for
alpar@9	1925 each column, and places all columns in the degree lists. Not user-callable.
alpar@9	1926 */
alpar@9	1927
alpar@9	1928 PRIVATE void init_scoring
alpar@9	1929 (
alpar@9	1930 /* === Parameters ======================================================= */
alpar@9	1931
alpar@9	1932 Int n_row, /* number of rows of A */
alpar@9	1933 Int n_col, /* number of columns of A */
alpar@9	1934 Colamd_Row Row [], /* of size n_row+1 */
alpar@9	1935 Colamd_Col Col [], /* of size n_col+1 */
alpar@9	1936 Int A [], /* column form and row form of A */
alpar@9	1937 Int head [], /* of size n_col+1 */
alpar@9	1938 double knobs [COLAMD_KNOBS],/* parameters */
alpar@9	1939 Int p_n_row2, / number of non-dense, non-empty rows */
alpar@9	1940 Int p_n_col2, / number of non-dense, non-empty columns */
alpar@9	1941 Int p_max_deg / maximum row degree */
alpar@9	1942 )
alpar@9	1943 {
alpar@9	1944 /* === Local variables ================================================== */
alpar@9	1945
alpar@9	1946 Int c ; /* a column index */
alpar@9	1947 Int r, row ; /* a row index */
alpar@9	1948 Int cp ; / a column pointer */
alpar@9	1949 Int deg ; /* degree of a row or column */
alpar@9	1950 Int cp_end ; / a pointer to the end of a column */
alpar@9	1951 Int new_cp ; / new column pointer */
alpar@9	1952 Int col_length ; /* length of pruned column */
alpar@9	1953 Int score ; /* current column score */
alpar@9	1954 Int n_col2 ; /* number of non-dense, non-empty columns */
alpar@9	1955 Int n_row2 ; /* number of non-dense, non-empty rows */
alpar@9	1956 Int dense_row_count ; /* remove rows with more entries than this */
alpar@9	1957 Int dense_col_count ; /* remove cols with more entries than this */
alpar@9	1958 Int min_score ; /* smallest column score */
alpar@9	1959 Int max_deg ; /* maximum row degree */
alpar@9	1960 Int next_col ; /* Used to add to degree list.*/
alpar@9	1961
alpar@9	1962 #ifndef NDEBUG
alpar@9	1963 Int debug_count ; /* debug only. */
alpar@9	1964 #endif /* NDEBUG */
alpar@9	1965
alpar@9	1966 /* === Extract knobs ==================================================== */
alpar@9	1967
alpar@9	1968 /* Note: if knobs contains a NaN, this is undefined: */
alpar@9	1969 if (knobs [COLAMD_DENSE_ROW] < 0)
alpar@9	1970 {
alpar@9	1971 /* only remove completely dense rows */
alpar@9	1972 dense_row_count = n_col-1 ;
alpar@9	1973 }
alpar@9	1974 else
alpar@9	1975 {
alpar@9	1976 dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
alpar@9	1977 }
alpar@9	1978 if (knobs [COLAMD_DENSE_COL] < 0)
alpar@9	1979 {
alpar@9	1980 /* only remove completely dense columns */
alpar@9	1981 dense_col_count = n_row-1 ;
alpar@9	1982 }
alpar@9	1983 else
alpar@9	1984 {
alpar@9	1985 dense_col_count =
alpar@9	1986 DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
alpar@9	1987 }
alpar@9	1988
alpar@9	1989 DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
alpar@9	1990 max_deg = 0 ;
alpar@9	1991 n_col2 = n_col ;
alpar@9	1992 n_row2 = n_row ;
alpar@9	1993
alpar@9	1994 /* === Kill empty columns =============================================== */
alpar@9	1995
alpar@9	1996 /* Put the empty columns at the end in their natural order, so that LU */
alpar@9	1997 /* factorization can proceed as far as possible. */
alpar@9	1998 for (c = n_col-1 ; c >= 0 ; c--)
alpar@9	1999 {
alpar@9	2000 deg = Col [c].length ;
alpar@9	2001 if (deg == 0)
alpar@9	2002 {
alpar@9	2003 /* this is a empty column, kill and order it last */
alpar@9	2004 Col [c].shared2.order = --n_col2 ;
alpar@9	2005 KILL_PRINCIPAL_COL (c) ;
alpar@9	2006 }
alpar@9	2007 }
alpar@9	2008 DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
alpar@9	2009
alpar@9	2010 /* === Kill dense columns =============================================== */
alpar@9	2011
alpar@9	2012 /* Put the dense columns at the end, in their natural order */
alpar@9	2013 for (c = n_col-1 ; c >= 0 ; c--)
alpar@9	2014 {
alpar@9	2015 /* skip any dead columns */
alpar@9	2016 if (COL_IS_DEAD (c))
alpar@9	2017 {
alpar@9	2018 continue ;
alpar@9	2019 }
alpar@9	2020 deg = Col [c].length ;
alpar@9	2021 if (deg > dense_col_count)
alpar@9	2022 {
alpar@9	2023 /* this is a dense column, kill and order it last */
alpar@9	2024 Col [c].shared2.order = --n_col2 ;
alpar@9	2025 /* decrement the row degrees */
alpar@9	2026 cp = &A [Col [c].start] ;
alpar@9	2027 cp_end = cp + Col [c].length ;
alpar@9	2028 while (cp < cp_end)
alpar@9	2029 {
alpar@9	2030 Row [*cp++].shared1.degree-- ;
alpar@9	2031 }
alpar@9	2032 KILL_PRINCIPAL_COL (c) ;
alpar@9	2033 }
alpar@9	2034 }
alpar@9	2035 DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
alpar@9	2036
alpar@9	2037 /* === Kill dense and empty rows ======================================== */
alpar@9	2038
alpar@9	2039 for (r = 0 ; r < n_row ; r++)
alpar@9	2040 {
alpar@9	2041 deg = Row [r].shared1.degree ;
alpar@9	2042 ASSERT (deg >= 0 && deg <= n_col) ;
alpar@9	2043 if (deg > dense_row_count \|\| deg == 0)
alpar@9	2044 {
alpar@9	2045 /* kill a dense or empty row */
alpar@9	2046 KILL_ROW (r) ;
alpar@9	2047 --n_row2 ;
alpar@9	2048 }
alpar@9	2049 else
alpar@9	2050 {
alpar@9	2051 /* keep track of max degree of remaining rows */
alpar@9	2052 max_deg = MAX (max_deg, deg) ;
alpar@9	2053 }
alpar@9	2054 }
alpar@9	2055 DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
alpar@9	2056
alpar@9	2057 /* === Compute initial column scores ==================================== */
alpar@9	2058
alpar@9	2059 /* At this point the row degrees are accurate. They reflect the number */
alpar@9	2060 /* of "live" (non-dense) columns in each row. No empty rows exist. */
alpar@9	2061 /* Some "live" columns may contain only dead rows, however. These are */
alpar@9	2062 /* pruned in the code below. */
alpar@9	2063
alpar@9	2064 /* now find the initial matlab score for each column */
alpar@9	2065 for (c = n_col-1 ; c >= 0 ; c--)
alpar@9	2066 {
alpar@9	2067 /* skip dead column */
alpar@9	2068 if (COL_IS_DEAD (c))
alpar@9	2069 {
alpar@9	2070 continue ;
alpar@9	2071 }
alpar@9	2072 score = 0 ;
alpar@9	2073 cp = &A [Col [c].start] ;
alpar@9	2074 new_cp = cp ;
alpar@9	2075 cp_end = cp + Col [c].length ;
alpar@9	2076 while (cp < cp_end)
alpar@9	2077 {
alpar@9	2078 /* get a row */
alpar@9	2079 row = *cp++ ;
alpar@9	2080 /* skip if dead */
alpar@9	2081 if (ROW_IS_DEAD (row))
alpar@9	2082 {
alpar@9	2083 continue ;
alpar@9	2084 }
alpar@9	2085 /* compact the column */
alpar@9	2086 *new_cp++ = row ;
alpar@9	2087 /* add row's external degree */
alpar@9	2088 score += Row [row].shared1.degree - 1 ;
alpar@9	2089 /* guard against integer overflow */
alpar@9	2090 score = MIN (score, n_col) ;
alpar@9	2091 }
alpar@9	2092 /* determine pruned column length */
alpar@9	2093 col_length = (Int) (new_cp - &A [Col [c].start]) ;
alpar@9	2094 if (col_length == 0)
alpar@9	2095 {
alpar@9	2096 /* a newly-made null column (all rows in this col are "dense" */
alpar@9	2097 /* and have already been killed) */
alpar@9	2098 DEBUG2 (("Newly null killed: %d\n", c)) ;
alpar@9	2099 Col [c].shared2.order = --n_col2 ;
alpar@9	2100 KILL_PRINCIPAL_COL (c) ;
alpar@9	2101 }
alpar@9	2102 else
alpar@9	2103 {
alpar@9	2104 /* set column length and set score */
alpar@9	2105 ASSERT (score >= 0) ;
alpar@9	2106 ASSERT (score <= n_col) ;
alpar@9	2107 Col [c].length = col_length ;
alpar@9	2108 Col [c].shared2.score = score ;
alpar@9	2109 }
alpar@9	2110 }
alpar@9	2111 DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
alpar@9	2112 n_col-n_col2)) ;
alpar@9	2113
alpar@9	2114 /* At this point, all empty rows and columns are dead. All live columns */
alpar@9	2115 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
alpar@9	2116 /* yet). Rows may contain dead columns, but all live rows contain at */
alpar@9	2117 /* least one live column. */
alpar@9	2118
alpar@9	2119 #ifndef NDEBUG
alpar@9	2120 debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
alpar@9	2121 #endif /* NDEBUG */
alpar@9	2122
alpar@9	2123 /* === Initialize degree lists ========================================== */
alpar@9	2124
alpar@9	2125 #ifndef NDEBUG
alpar@9	2126 debug_count = 0 ;
alpar@9	2127 #endif /* NDEBUG */
alpar@9	2128
alpar@9	2129 /* clear the hash buckets */
alpar@9	2130 for (c = 0 ; c <= n_col ; c++)
alpar@9	2131 {
alpar@9	2132 head [c] = EMPTY ;
alpar@9	2133 }
alpar@9	2134 min_score = n_col ;
alpar@9	2135 /* place in reverse order, so low column indices are at the front */
alpar@9	2136 /* of the lists. This is to encourage natural tie-breaking */
alpar@9	2137 for (c = n_col-1 ; c >= 0 ; c--)
alpar@9	2138 {
alpar@9	2139 /* only add principal columns to degree lists */
alpar@9	2140 if (COL_IS_ALIVE (c))
alpar@9	2141 {
alpar@9	2142 DEBUG4 (("place %d score %d minscore %d ncol %d\n",
alpar@9	2143 c, Col [c].shared2.score, min_score, n_col)) ;
alpar@9	2144
alpar@9	2145 /* === Add columns score to DList =============================== */
alpar@9	2146
alpar@9	2147 score = Col [c].shared2.score ;
alpar@9	2148
alpar@9	2149 ASSERT (min_score >= 0) ;
alpar@9	2150 ASSERT (min_score <= n_col) ;
alpar@9	2151 ASSERT (score >= 0) ;
alpar@9	2152 ASSERT (score <= n_col) ;
alpar@9	2153 ASSERT (head [score] >= EMPTY) ;
alpar@9	2154
alpar@9	2155 /* now add this column to dList at proper score location */
alpar@9	2156 next_col = head [score] ;
alpar@9	2157 Col [c].shared3.prev = EMPTY ;
alpar@9	2158 Col [c].shared4.degree_next = next_col ;
alpar@9	2159
alpar@9	2160 /* if there already was a column with the same score, set its */
alpar@9	2161 /* previous pointer to this new column */
alpar@9	2162 if (next_col != EMPTY)
alpar@9	2163 {
alpar@9	2164 Col [next_col].shared3.prev = c ;
alpar@9	2165 }
alpar@9	2166 head [score] = c ;
alpar@9	2167
alpar@9	2168 /* see if this score is less than current min */
alpar@9	2169 min_score = MIN (min_score, score) ;
alpar@9	2170
alpar@9	2171 #ifndef NDEBUG
alpar@9	2172 debug_count++ ;
alpar@9	2173 #endif /* NDEBUG */
alpar@9	2174
alpar@9	2175 }
alpar@9	2176 }
alpar@9	2177
alpar@9	2178 #ifndef NDEBUG
alpar@9	2179 DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
alpar@9	2180 debug_count, n_col, n_col-debug_count)) ;
alpar@9	2181 ASSERT (debug_count == n_col2) ;
alpar@9	2182 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
alpar@9	2183 #endif /* NDEBUG */
alpar@9	2184
alpar@9	2185 /* === Return number of remaining columns, and max row degree =========== */
alpar@9	2186
alpar@9	2187 *p_n_col2 = n_col2 ;
alpar@9	2188 *p_n_row2 = n_row2 ;
alpar@9	2189 *p_max_deg = max_deg ;
alpar@9	2190 }
alpar@9	2191
alpar@9	2192
alpar@9	2193 /* ========================================================================== */
alpar@9	2194 /* === find_ordering ======================================================== */
alpar@9	2195 /* ========================================================================== */
alpar@9	2196
alpar@9	2197 /*
alpar@9	2198 Order the principal columns of the supercolumn form of the matrix
alpar@9	2199 (no supercolumns on input). Uses a minimum approximate column minimum
alpar@9	2200 degree ordering method. Not user-callable.
alpar@9	2201 */
alpar@9	2202
alpar@9	2203 PRIVATE Int find_ordering /* return the number of garbage collections */
alpar@9	2204 (
alpar@9	2205 /* === Parameters ======================================================= */
alpar@9	2206
alpar@9	2207 Int n_row, /* number of rows of A */
alpar@9	2208 Int n_col, /* number of columns of A */
alpar@9	2209 Int Alen, /* size of A, 2nnz + n_col or larger /
alpar@9	2210 Colamd_Row Row [], /* of size n_row+1 */
alpar@9	2211 Colamd_Col Col [], /* of size n_col+1 */
alpar@9	2212 Int A [], /* column form and row form of A */
alpar@9	2213 Int head [], /* of size n_col+1 */
alpar@9	2214 Int n_col2, /* Remaining columns to order */
alpar@9	2215 Int max_deg, /* Maximum row degree */
alpar@9	2216 Int pfree, /* index of first free slot (2nnz on entry) /
alpar@9	2217 Int aggressive
alpar@9	2218 )
alpar@9	2219 {
alpar@9	2220 /* === Local variables ================================================== */
alpar@9	2221
alpar@9	2222 Int k ; /* current pivot ordering step */
alpar@9	2223 Int pivot_col ; /* current pivot column */
alpar@9	2224 Int cp ; / a column pointer */
alpar@9	2225 Int rp ; / a row pointer */
alpar@9	2226 Int pivot_row ; /* current pivot row */
alpar@9	2227 Int new_cp ; / modified column pointer */
alpar@9	2228 Int new_rp ; / modified row pointer */
alpar@9	2229 Int pivot_row_start ; /* pointer to start of pivot row */
alpar@9	2230 Int pivot_row_degree ; /* number of columns in pivot row */
alpar@9	2231 Int pivot_row_length ; /* number of supercolumns in pivot row */
alpar@9	2232 Int pivot_col_score ; /* score of pivot column */
alpar@9	2233 Int needed_memory ; /* free space needed for pivot row */
alpar@9	2234 Int cp_end ; / pointer to the end of a column */
alpar@9	2235 Int rp_end ; / pointer to the end of a row */
alpar@9	2236 Int row ; /* a row index */
alpar@9	2237 Int col ; /* a column index */
alpar@9	2238 Int max_score ; /* maximum possible score */
alpar@9	2239 Int cur_score ; /* score of current column */
alpar@9	2240 unsigned Int hash ; /* hash value for supernode detection */
alpar@9	2241 Int head_column ; /* head of hash bucket */
alpar@9	2242 Int first_col ; /* first column in hash bucket */
alpar@9	2243 Int tag_mark ; /* marker value for mark array */
alpar@9	2244 Int row_mark ; /* Row [row].shared2.mark */
alpar@9	2245 Int set_difference ; /* set difference size of row with pivot row */
alpar@9	2246 Int min_score ; /* smallest column score */
alpar@9	2247 Int col_thickness ; /* "thickness" (no. of columns in a supercol) */
alpar@9	2248 Int max_mark ; /* maximum value of tag_mark */
alpar@9	2249 Int pivot_col_thickness ; /* number of columns represented by pivot col */
alpar@9	2250 Int prev_col ; /* Used by Dlist operations. */
alpar@9	2251 Int next_col ; /* Used by Dlist operations. */
alpar@9	2252 Int ngarbage ; /* number of garbage collections performed */
alpar@9	2253
alpar@9	2254 #ifndef NDEBUG
alpar@9	2255 Int debug_d ; /* debug loop counter */
alpar@9	2256 Int debug_step = 0 ; /* debug loop counter */
alpar@9	2257 #endif /* NDEBUG */
alpar@9	2258
alpar@9	2259 /* === Initialization and clear mark ==================================== */
alpar@9	2260
alpar@9	2261 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */
alpar@9	2262 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
alpar@9	2263 min_score = 0 ;
alpar@9	2264 ngarbage = 0 ;
alpar@9	2265 DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
alpar@9	2266
alpar@9	2267 /* === Order the columns ================================================ */
alpar@9	2268
alpar@9	2269 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
alpar@9	2270 {
alpar@9	2271
alpar@9	2272 #ifndef NDEBUG
alpar@9	2273 if (debug_step % 100 == 0)
alpar@9	2274 {
alpar@9	2275 DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ;
alpar@9	2276 }
alpar@9	2277 else
alpar@9	2278 {
alpar@9	2279 DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
alpar@9	2280 }
alpar@9	2281 debug_step++ ;
alpar@9	2282 debug_deg_lists (n_row, n_col, Row, Col, head,
alpar@9	2283 min_score, n_col2-k, max_deg) ;
alpar@9	2284 debug_matrix (n_row, n_col, Row, Col, A) ;
alpar@9	2285 #endif /* NDEBUG */
alpar@9	2286
alpar@9	2287 /* === Select pivot column, and order it ============================ */
alpar@9	2288
alpar@9	2289 /* make sure degree list isn't empty */
alpar@9	2290 ASSERT (min_score >= 0) ;
alpar@9	2291 ASSERT (min_score <= n_col) ;
alpar@9	2292 ASSERT (head [min_score] >= EMPTY) ;
alpar@9	2293
alpar@9	2294 #ifndef NDEBUG
alpar@9	2295 for (debug_d = 0 ; debug_d < min_score ; debug_d++)
alpar@9	2296 {
alpar@9	2297 ASSERT (head [debug_d] == EMPTY) ;
alpar@9	2298 }
alpar@9	2299 #endif /* NDEBUG */
alpar@9	2300
alpar@9	2301 /* get pivot column from head of minimum degree list */
alpar@9	2302 while (head [min_score] == EMPTY && min_score < n_col)
alpar@9	2303 {
alpar@9	2304 min_score++ ;
alpar@9	2305 }
alpar@9	2306 pivot_col = head [min_score] ;
alpar@9	2307 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
alpar@9	2308 next_col = Col [pivot_col].shared4.degree_next ;
alpar@9	2309 head [min_score] = next_col ;
alpar@9	2310 if (next_col != EMPTY)
alpar@9	2311 {
alpar@9	2312 Col [next_col].shared3.prev = EMPTY ;
alpar@9	2313 }
alpar@9	2314
alpar@9	2315 ASSERT (COL_IS_ALIVE (pivot_col)) ;
alpar@9	2316
alpar@9	2317 /* remember score for defrag check */
alpar@9	2318 pivot_col_score = Col [pivot_col].shared2.score ;
alpar@9	2319
alpar@9	2320 /* the pivot column is the kth column in the pivot order */
alpar@9	2321 Col [pivot_col].shared2.order = k ;
alpar@9	2322
alpar@9	2323 /* increment order count by column thickness */
alpar@9	2324 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
alpar@9	2325 k += pivot_col_thickness ;
alpar@9	2326 ASSERT (pivot_col_thickness > 0) ;
alpar@9	2327 DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
alpar@9	2328
alpar@9	2329 /* === Garbage_collection, if necessary ============================= */
alpar@9	2330
alpar@9	2331 needed_memory = MIN (pivot_col_score, n_col - k) ;
alpar@9	2332 if (pfree + needed_memory >= Alen)
alpar@9	2333 {
alpar@9	2334 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
alpar@9	2335 ngarbage++ ;
alpar@9	2336 /* after garbage collection we will have enough */
alpar@9	2337 ASSERT (pfree + needed_memory < Alen) ;
alpar@9	2338 /* garbage collection has wiped out the Row[].shared2.mark array */
alpar@9	2339 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
alpar@9	2340
alpar@9	2341 #ifndef NDEBUG
alpar@9	2342 debug_matrix (n_row, n_col, Row, Col, A) ;
alpar@9	2343 #endif /* NDEBUG */
alpar@9	2344 }
alpar@9	2345
alpar@9	2346 /* === Compute pivot row pattern ==================================== */
alpar@9	2347
alpar@9	2348 /* get starting location for this new merged row */
alpar@9	2349 pivot_row_start = pfree ;
alpar@9	2350
alpar@9	2351 /* initialize new row counts to zero */
alpar@9	2352 pivot_row_degree = 0 ;
alpar@9	2353
alpar@9	2354 /* tag pivot column as having been visited so it isn't included */
alpar@9	2355 /* in merged pivot row */
alpar@9	2356 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
alpar@9	2357
alpar@9	2358 /* pivot row is the union of all rows in the pivot column pattern */
alpar@9	2359 cp = &A [Col [pivot_col].start] ;
alpar@9	2360 cp_end = cp + Col [pivot_col].length ;
alpar@9	2361 while (cp < cp_end)
alpar@9	2362 {
alpar@9	2363 /* get a row */
alpar@9	2364 row = *cp++ ;
alpar@9	2365 DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
alpar@9	2366 /* skip if row is dead */
alpar@9	2367 if (ROW_IS_ALIVE (row))
alpar@9	2368 {
alpar@9	2369 rp = &A [Row [row].start] ;
alpar@9	2370 rp_end = rp + Row [row].length ;
alpar@9	2371 while (rp < rp_end)
alpar@9	2372 {
alpar@9	2373 /* get a column */
alpar@9	2374 col = *rp++ ;
alpar@9	2375 /* add the column, if alive and untagged */
alpar@9	2376 col_thickness = Col [col].shared1.thickness ;
alpar@9	2377 if (col_thickness > 0 && COL_IS_ALIVE (col))
alpar@9	2378 {
alpar@9	2379 /* tag column in pivot row */
alpar@9	2380 Col [col].shared1.thickness = -col_thickness ;
alpar@9	2381 ASSERT (pfree < Alen) ;
alpar@9	2382 /* place column in pivot row */
alpar@9	2383 A [pfree++] = col ;
alpar@9	2384 pivot_row_degree += col_thickness ;
alpar@9	2385 }
alpar@9	2386 }
alpar@9	2387 }
alpar@9	2388 }
alpar@9	2389
alpar@9	2390 /* clear tag on pivot column */
alpar@9	2391 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
alpar@9	2392 max_deg = MAX (max_deg, pivot_row_degree) ;
alpar@9	2393
alpar@9	2394 #ifndef NDEBUG
alpar@9	2395 DEBUG3 (("check2\n")) ;
alpar@9	2396 debug_mark (n_row, Row, tag_mark, max_mark) ;
alpar@9	2397 #endif /* NDEBUG */
alpar@9	2398
alpar@9	2399 /* === Kill all rows used to construct pivot row ==================== */
alpar@9	2400
alpar@9	2401 /* also kill pivot row, temporarily */
alpar@9	2402 cp = &A [Col [pivot_col].start] ;
alpar@9	2403 cp_end = cp + Col [pivot_col].length ;
alpar@9	2404 while (cp < cp_end)
alpar@9	2405 {
alpar@9	2406 /* may be killing an already dead row */
alpar@9	2407 row = *cp++ ;
alpar@9	2408 DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
alpar@9	2409 KILL_ROW (row) ;
alpar@9	2410 }
alpar@9	2411
alpar@9	2412 /* === Select a row index to use as the new pivot row =============== */
alpar@9	2413
alpar@9	2414 pivot_row_length = pfree - pivot_row_start ;
alpar@9	2415 if (pivot_row_length > 0)
alpar@9	2416 {
alpar@9	2417 /* pick the "pivot" row arbitrarily (first row in col) */
alpar@9	2418 pivot_row = A [Col [pivot_col].start] ;
alpar@9	2419 DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
alpar@9	2420 }
alpar@9	2421 else
alpar@9	2422 {
alpar@9	2423 /* there is no pivot row, since it is of zero length */
alpar@9	2424 pivot_row = EMPTY ;
alpar@9	2425 ASSERT (pivot_row_length == 0) ;
alpar@9	2426 }
alpar@9	2427 ASSERT (Col [pivot_col].length > 0 \|\| pivot_row_length == 0) ;
alpar@9	2428
alpar@9	2429 /* === Approximate degree computation =============================== */
alpar@9	2430
alpar@9	2431 /* Here begins the computation of the approximate degree. The column */
alpar@9	2432 /* score is the sum of the pivot row "length", plus the size of the */
alpar@9	2433 /* set differences of each row in the column minus the pattern of the */
alpar@9	2434 /* pivot row itself. The column ("thickness") itself is also */
alpar@9	2435 /* excluded from the column score (we thus use an approximate */
alpar@9	2436 /* external degree). */
alpar@9	2437
alpar@9	2438 /* The time taken by the following code (compute set differences, and */
alpar@9	2439 /* add them up) is proportional to the size of the data structure */
alpar@9	2440 /* being scanned - that is, the sum of the sizes of each column in */
alpar@9	2441 /* the pivot row. Thus, the amortized time to compute a column score */
alpar@9	2442 /* is proportional to the size of that column (where size, in this */
alpar@9	2443 /* context, is the column "length", or the number of row indices */
alpar@9	2444 /* in that column). The number of row indices in a column is */
alpar@9	2445 /* monotonically non-decreasing, from the length of the original */
alpar@9	2446 /* column on input to colamd. */
alpar@9	2447
alpar@9	2448 /* === Compute set differences ====================================== */
alpar@9	2449
alpar@9	2450 DEBUG3 ((" Computing set differences phase. \n")) ;
alpar@9	2451
alpar@9	2452 /* pivot row is currently dead - it will be revived later. */
alpar@9	2453
alpar@9	2454 DEBUG3 (("Pivot row: ")) ;
alpar@9	2455 /* for each column in pivot row */
alpar@9	2456 rp = &A [pivot_row_start] ;
alpar@9	2457 rp_end = rp + pivot_row_length ;
alpar@9	2458 while (rp < rp_end)
alpar@9	2459 {
alpar@9	2460 col = *rp++ ;
alpar@9	2461 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
alpar@9	2462 DEBUG3 (("Col: %d\n", col)) ;
alpar@9	2463
alpar@9	2464 /* clear tags used to construct pivot row pattern */
alpar@9	2465 col_thickness = -Col [col].shared1.thickness ;
alpar@9	2466 ASSERT (col_thickness > 0) ;
alpar@9	2467 Col [col].shared1.thickness = col_thickness ;
alpar@9	2468
alpar@9	2469 /* === Remove column from degree list =========================== */
alpar@9	2470
alpar@9	2471 cur_score = Col [col].shared2.score ;
alpar@9	2472 prev_col = Col [col].shared3.prev ;
alpar@9	2473 next_col = Col [col].shared4.degree_next ;
alpar@9	2474 ASSERT (cur_score >= 0) ;
alpar@9	2475 ASSERT (cur_score <= n_col) ;
alpar@9	2476 ASSERT (cur_score >= EMPTY) ;
alpar@9	2477 if (prev_col == EMPTY)
alpar@9	2478 {
alpar@9	2479 head [cur_score] = next_col ;
alpar@9	2480 }
alpar@9	2481 else
alpar@9	2482 {
alpar@9	2483 Col [prev_col].shared4.degree_next = next_col ;
alpar@9	2484 }
alpar@9	2485 if (next_col != EMPTY)
alpar@9	2486 {
alpar@9	2487 Col [next_col].shared3.prev = prev_col ;
alpar@9	2488 }
alpar@9	2489
alpar@9	2490 /* === Scan the column ========================================== */
alpar@9	2491
alpar@9	2492 cp = &A [Col [col].start] ;
alpar@9	2493 cp_end = cp + Col [col].length ;
alpar@9	2494 while (cp < cp_end)
alpar@9	2495 {
alpar@9	2496 /* get a row */
alpar@9	2497 row = *cp++ ;
alpar@9	2498 row_mark = Row [row].shared2.mark ;
alpar@9	2499 /* skip if dead */
alpar@9	2500 if (ROW_IS_MARKED_DEAD (row_mark))
alpar@9	2501 {
alpar@9	2502 continue ;
alpar@9	2503 }
alpar@9	2504 ASSERT (row != pivot_row) ;
alpar@9	2505 set_difference = row_mark - tag_mark ;
alpar@9	2506 /* check if the row has been seen yet */
alpar@9	2507 if (set_difference < 0)
alpar@9	2508 {
alpar@9	2509 ASSERT (Row [row].shared1.degree <= max_deg) ;
alpar@9	2510 set_difference = Row [row].shared1.degree ;
alpar@9	2511 }
alpar@9	2512 /* subtract column thickness from this row's set difference */
alpar@9	2513 set_difference -= col_thickness ;
alpar@9	2514 ASSERT (set_difference >= 0) ;
alpar@9	2515 /* absorb this row if the set difference becomes zero */
alpar@9	2516 if (set_difference == 0 && aggressive)
alpar@9	2517 {
alpar@9	2518 DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
alpar@9	2519 KILL_ROW (row) ;
alpar@9	2520 }
alpar@9	2521 else
alpar@9	2522 {
alpar@9	2523 /* save the new mark */
alpar@9	2524 Row [row].shared2.mark = set_difference + tag_mark ;
alpar@9	2525 }
alpar@9	2526 }
alpar@9	2527 }
alpar@9	2528
alpar@9	2529 #ifndef NDEBUG
alpar@9	2530 debug_deg_lists (n_row, n_col, Row, Col, head,
alpar@9	2531 min_score, n_col2-k-pivot_row_degree, max_deg) ;
alpar@9	2532 #endif /* NDEBUG */
alpar@9	2533
alpar@9	2534 /* === Add up set differences for each column ======================= */
alpar@9	2535
alpar@9	2536 DEBUG3 ((" Adding set differences phase. \n")) ;
alpar@9	2537
alpar@9	2538 /* for each column in pivot row */
alpar@9	2539 rp = &A [pivot_row_start] ;
alpar@9	2540 rp_end = rp + pivot_row_length ;
alpar@9	2541 while (rp < rp_end)
alpar@9	2542 {
alpar@9	2543 /* get a column */
alpar@9	2544 col = *rp++ ;
alpar@9	2545 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
alpar@9	2546 hash = 0 ;
alpar@9	2547 cur_score = 0 ;
alpar@9	2548 cp = &A [Col [col].start] ;
alpar@9	2549 /* compact the column */
alpar@9	2550 new_cp = cp ;
alpar@9	2551 cp_end = cp + Col [col].length ;
alpar@9	2552
alpar@9	2553 DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
alpar@9	2554
alpar@9	2555 while (cp < cp_end)
alpar@9	2556 {
alpar@9	2557 /* get a row */
alpar@9	2558 row = *cp++ ;
alpar@9	2559 ASSERT(row >= 0 && row < n_row) ;
alpar@9	2560 row_mark = Row [row].shared2.mark ;
alpar@9	2561 /* skip if dead */
alpar@9	2562 if (ROW_IS_MARKED_DEAD (row_mark))
alpar@9	2563 {
alpar@9	2564 DEBUG4 ((" Row %d, dead\n", row)) ;
alpar@9	2565 continue ;
alpar@9	2566 }
alpar@9	2567 DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
alpar@9	2568 ASSERT (row_mark >= tag_mark) ;
alpar@9	2569 /* compact the column */
alpar@9	2570 *new_cp++ = row ;
alpar@9	2571 /* compute hash function */
alpar@9	2572 hash += row ;
alpar@9	2573 /* add set difference */
alpar@9	2574 cur_score += row_mark - tag_mark ;
alpar@9	2575 /* integer overflow... */
alpar@9	2576 cur_score = MIN (cur_score, n_col) ;
alpar@9	2577 }
alpar@9	2578
alpar@9	2579 /* recompute the column's length */
alpar@9	2580 Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
alpar@9	2581
alpar@9	2582 /* === Further mass elimination ================================= */
alpar@9	2583
alpar@9	2584 if (Col [col].length == 0)
alpar@9	2585 {
alpar@9	2586 DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
alpar@9	2587 /* nothing left but the pivot row in this column */
alpar@9	2588 KILL_PRINCIPAL_COL (col) ;
alpar@9	2589 pivot_row_degree -= Col [col].shared1.thickness ;
alpar@9	2590 ASSERT (pivot_row_degree >= 0) ;
alpar@9	2591 /* order it */
alpar@9	2592 Col [col].shared2.order = k ;
alpar@9	2593 /* increment order count by column thickness */
alpar@9	2594 k += Col [col].shared1.thickness ;
alpar@9	2595 }
alpar@9	2596 else
alpar@9	2597 {
alpar@9	2598 /* === Prepare for supercolumn detection ==================== */
alpar@9	2599
alpar@9	2600 DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
alpar@9	2601
alpar@9	2602 /* save score so far */
alpar@9	2603 Col [col].shared2.score = cur_score ;
alpar@9	2604
alpar@9	2605 /* add column to hash table, for supercolumn detection */
alpar@9	2606 hash %= n_col + 1 ;
alpar@9	2607
alpar@9	2608 DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
alpar@9	2609 ASSERT (((Int) hash) <= n_col) ;
alpar@9	2610
alpar@9	2611 head_column = head [hash] ;
alpar@9	2612 if (head_column > EMPTY)
alpar@9	2613 {
alpar@9	2614 /* degree list "hash" is non-empty, use prev (shared3) of */
alpar@9	2615 /* first column in degree list as head of hash bucket */
alpar@9	2616 first_col = Col [head_column].shared3.headhash ;
alpar@9	2617 Col [head_column].shared3.headhash = col ;
alpar@9	2618 }
alpar@9	2619 else
alpar@9	2620 {
alpar@9	2621 /* degree list "hash" is empty, use head as hash bucket */
alpar@9	2622 first_col = - (head_column + 2) ;
alpar@9	2623 head [hash] = - (col + 2) ;
alpar@9	2624 }
alpar@9	2625 Col [col].shared4.hash_next = first_col ;
alpar@9	2626
alpar@9	2627 /* save hash function in Col [col].shared3.hash */
alpar@9	2628 Col [col].shared3.hash = (Int) hash ;
alpar@9	2629 ASSERT (COL_IS_ALIVE (col)) ;
alpar@9	2630 }
alpar@9	2631 }
alpar@9	2632
alpar@9	2633 /* The approximate external column degree is now computed. */
alpar@9	2634
alpar@9	2635 /* === Supercolumn detection ======================================== */
alpar@9	2636
alpar@9	2637 DEBUG3 ((" Supercolumn detection phase. \n")) ;
alpar@9	2638
alpar@9	2639 detect_super_cols (
alpar@9	2640
alpar@9	2641 #ifndef NDEBUG
alpar@9	2642 n_col, Row,
alpar@9	2643 #endif /* NDEBUG */
alpar@9	2644
alpar@9	2645 Col, A, head, pivot_row_start, pivot_row_length) ;
alpar@9	2646
alpar@9	2647 /* === Kill the pivotal column ====================================== */
alpar@9	2648
alpar@9	2649 KILL_PRINCIPAL_COL (pivot_col) ;
alpar@9	2650
alpar@9	2651 /* === Clear mark =================================================== */
alpar@9	2652
alpar@9	2653 tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
alpar@9	2654
alpar@9	2655 #ifndef NDEBUG
alpar@9	2656 DEBUG3 (("check3\n")) ;
alpar@9	2657 debug_mark (n_row, Row, tag_mark, max_mark) ;
alpar@9	2658 #endif /* NDEBUG */
alpar@9	2659
alpar@9	2660 /* === Finalize the new pivot row, and column scores ================ */
alpar@9	2661
alpar@9	2662 DEBUG3 ((" Finalize scores phase. \n")) ;
alpar@9	2663
alpar@9	2664 /* for each column in pivot row */
alpar@9	2665 rp = &A [pivot_row_start] ;
alpar@9	2666 /* compact the pivot row */
alpar@9	2667 new_rp = rp ;
alpar@9	2668 rp_end = rp + pivot_row_length ;
alpar@9	2669 while (rp < rp_end)
alpar@9	2670 {
alpar@9	2671 col = *rp++ ;
alpar@9	2672 /* skip dead columns */
alpar@9	2673 if (COL_IS_DEAD (col))
alpar@9	2674 {
alpar@9	2675 continue ;
alpar@9	2676 }
alpar@9	2677 *new_rp++ = col ;
alpar@9	2678 /* add new pivot row to column */
alpar@9	2679 A [Col [col].start + (Col [col].length++)] = pivot_row ;
alpar@9	2680
alpar@9	2681 /* retrieve score so far and add on pivot row's degree. */
alpar@9	2682 /* (we wait until here for this in case the pivot */
alpar@9	2683 /* row's degree was reduced due to mass elimination). */
alpar@9	2684 cur_score = Col [col].shared2.score + pivot_row_degree ;
alpar@9	2685
alpar@9	2686 /* calculate the max possible score as the number of */
alpar@9	2687 /* external columns minus the 'k' value minus the */
alpar@9	2688 /* columns thickness */
alpar@9	2689 max_score = n_col - k - Col [col].shared1.thickness ;
alpar@9	2690
alpar@9	2691 /* make the score the external degree of the union-of-rows */
alpar@9	2692 cur_score -= Col [col].shared1.thickness ;
alpar@9	2693
alpar@9	2694 /* make sure score is less or equal than the max score */
alpar@9	2695 cur_score = MIN (cur_score, max_score) ;
alpar@9	2696 ASSERT (cur_score >= 0) ;
alpar@9	2697
alpar@9	2698 /* store updated score */
alpar@9	2699 Col [col].shared2.score = cur_score ;
alpar@9	2700
alpar@9	2701 /* === Place column back in degree list ========================= */
alpar@9	2702
alpar@9	2703 ASSERT (min_score >= 0) ;
alpar@9	2704 ASSERT (min_score <= n_col) ;
alpar@9	2705 ASSERT (cur_score >= 0) ;
alpar@9	2706 ASSERT (cur_score <= n_col) ;
alpar@9	2707 ASSERT (head [cur_score] >= EMPTY) ;
alpar@9	2708 next_col = head [cur_score] ;
alpar@9	2709 Col [col].shared4.degree_next = next_col ;
alpar@9	2710 Col [col].shared3.prev = EMPTY ;
alpar@9	2711 if (next_col != EMPTY)
alpar@9	2712 {
alpar@9	2713 Col [next_col].shared3.prev = col ;
alpar@9	2714 }
alpar@9	2715 head [cur_score] = col ;
alpar@9	2716
alpar@9	2717 /* see if this score is less than current min */
alpar@9	2718 min_score = MIN (min_score, cur_score) ;
alpar@9	2719
alpar@9	2720 }
alpar@9	2721
alpar@9	2722 #ifndef NDEBUG
alpar@9	2723 debug_deg_lists (n_row, n_col, Row, Col, head,
alpar@9	2724 min_score, n_col2-k, max_deg) ;
alpar@9	2725 #endif /* NDEBUG */
alpar@9	2726
alpar@9	2727 /* === Resurrect the new pivot row ================================== */
alpar@9	2728
alpar@9	2729 if (pivot_row_degree > 0)
alpar@9	2730 {
alpar@9	2731 /* update pivot row length to reflect any cols that were killed */
alpar@9	2732 /* during super-col detection and mass elimination */
alpar@9	2733 Row [pivot_row].start = pivot_row_start ;
alpar@9	2734 Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
alpar@9	2735 ASSERT (Row [pivot_row].length > 0) ;
alpar@9	2736 Row [pivot_row].shared1.degree = pivot_row_degree ;
alpar@9	2737 Row [pivot_row].shared2.mark = 0 ;
alpar@9	2738 /* pivot row is no longer dead */
alpar@9	2739
alpar@9	2740 DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
alpar@9	2741 pivot_row, pivot_row_degree)) ;
alpar@9	2742 }
alpar@9	2743 }
alpar@9	2744
alpar@9	2745 /* === All principal columns have now been ordered ====================== */
alpar@9	2746
alpar@9	2747 return (ngarbage) ;
alpar@9	2748 }
alpar@9	2749
alpar@9	2750
alpar@9	2751 /* ========================================================================== */
alpar@9	2752 /* === order_children ======================================================= */
alpar@9	2753 /* ========================================================================== */
alpar@9	2754
alpar@9	2755 /*
alpar@9	2756 The find_ordering routine has ordered all of the principal columns (the
alpar@9	2757 representatives of the supercolumns). The non-principal columns have not
alpar@9	2758 yet been ordered. This routine orders those columns by walking up the
alpar@9	2759 parent tree (a column is a child of the column which absorbed it). The
alpar@9	2760 final permutation vector is then placed in p [0 ... n_col-1], with p [0]
alpar@9	2761 being the first column, and p [n_col-1] being the last. It doesn't look
alpar@9	2762 like it at first glance, but be assured that this routine takes time linear
alpar@9	2763 in the number of columns. Although not immediately obvious, the time
alpar@9	2764 taken by this routine is O (n_col), that is, linear in the number of
alpar@9	2765 columns. Not user-callable.
alpar@9	2766 */
alpar@9	2767
alpar@9	2768 PRIVATE void order_children
alpar@9	2769 (
alpar@9	2770 /* === Parameters ======================================================= */
alpar@9	2771
alpar@9	2772 Int n_col, /* number of columns of A */
alpar@9	2773 Colamd_Col Col [], /* of size n_col+1 */
alpar@9	2774 Int p [] /* p [0 ... n_col-1] is the column permutation*/
alpar@9	2775 )
alpar@9	2776 {
alpar@9	2777 /* === Local variables ================================================== */
alpar@9	2778
alpar@9	2779 Int i ; /* loop counter for all columns */
alpar@9	2780 Int c ; /* column index */
alpar@9	2781 Int parent ; /* index of column's parent */
alpar@9	2782 Int order ; /* column's order */
alpar@9	2783
alpar@9	2784 /* === Order each non-principal column ================================== */
alpar@9	2785
alpar@9	2786 for (i = 0 ; i < n_col ; i++)
alpar@9	2787 {
alpar@9	2788 /* find an un-ordered non-principal column */
alpar@9	2789 ASSERT (COL_IS_DEAD (i)) ;
alpar@9	2790 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
alpar@9	2791 {
alpar@9	2792 parent = i ;
alpar@9	2793 /* once found, find its principal parent */
alpar@9	2794 do
alpar@9	2795 {
alpar@9	2796 parent = Col [parent].shared1.parent ;
alpar@9	2797 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
alpar@9	2798
alpar@9	2799 /* now, order all un-ordered non-principal columns along path */
alpar@9	2800 /* to this parent. collapse tree at the same time */
alpar@9	2801 c = i ;
alpar@9	2802 /* get order of parent */
alpar@9	2803 order = Col [parent].shared2.order ;
alpar@9	2804
alpar@9	2805 do
alpar@9	2806 {
alpar@9	2807 ASSERT (Col [c].shared2.order == EMPTY) ;
alpar@9	2808
alpar@9	2809 /* order this column */
alpar@9	2810 Col [c].shared2.order = order++ ;
alpar@9	2811 /* collaps tree */
alpar@9	2812 Col [c].shared1.parent = parent ;
alpar@9	2813
alpar@9	2814 /* get immediate parent of this column */
alpar@9	2815 c = Col [c].shared1.parent ;
alpar@9	2816
alpar@9	2817 /* continue until we hit an ordered column. There are */
alpar@9	2818 /* guarranteed not to be anymore unordered columns */
alpar@9	2819 /* above an ordered column */
alpar@9	2820 } while (Col [c].shared2.order == EMPTY) ;
alpar@9	2821
alpar@9	2822 /* re-order the super_col parent to largest order for this group */
alpar@9	2823 Col [parent].shared2.order = order ;
alpar@9	2824 }
alpar@9	2825 }
alpar@9	2826
alpar@9	2827 /* === Generate the permutation ========================================= */
alpar@9	2828
alpar@9	2829 for (c = 0 ; c < n_col ; c++)
alpar@9	2830 {
alpar@9	2831 p [Col [c].shared2.order] = c ;
alpar@9	2832 }
alpar@9	2833 }
alpar@9	2834
alpar@9	2835
alpar@9	2836 /* ========================================================================== */
alpar@9	2837 /* === detect_super_cols ==================================================== */
alpar@9	2838 /* ========================================================================== */
alpar@9	2839
alpar@9	2840 /*
alpar@9	2841 Detects supercolumns by finding matches between columns in the hash buckets.
alpar@9	2842 Check amongst columns in the set A [row_start ... row_start + row_length-1].
alpar@9	2843 The columns under consideration are currently not in the degree lists,
alpar@9	2844 and have already been placed in the hash buckets.
alpar@9	2845
alpar@9	2846 The hash bucket for columns whose hash function is equal to h is stored
alpar@9	2847 as follows:
alpar@9	2848
alpar@9	2849 if head [h] is >= 0, then head [h] contains a degree list, so:
alpar@9	2850
alpar@9	2851 head [h] is the first column in degree bucket h.
alpar@9	2852 Col [head [h]].headhash gives the first column in hash bucket h.
alpar@9	2853
alpar@9	2854 otherwise, the degree list is empty, and:
alpar@9	2855
alpar@9	2856 -(head [h] + 2) is the first column in hash bucket h.
alpar@9	2857
alpar@9	2858 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
alpar@9	2859 column" pointer. Col [c].shared3.hash is used instead as the hash number
alpar@9	2860 for that column. The value of Col [c].shared4.hash_next is the next column
alpar@9	2861 in the same hash bucket.
alpar@9	2862
alpar@9	2863 Assuming no, or "few" hash collisions, the time taken by this routine is
alpar@9	2864 linear in the sum of the sizes (lengths) of each column whose score has
alpar@9	2865 just been computed in the approximate degree computation.
alpar@9	2866 Not user-callable.
alpar@9	2867 */
alpar@9	2868
alpar@9	2869 PRIVATE void detect_super_cols
alpar@9	2870 (
alpar@9	2871 /* === Parameters ======================================================= */
alpar@9	2872
alpar@9	2873 #ifndef NDEBUG
alpar@9	2874 /* these two parameters are only needed when debugging is enabled: */
alpar@9	2875 Int n_col, /* number of columns of A */
alpar@9	2876 Colamd_Row Row [], /* of size n_row+1 */
alpar@9	2877 #endif /* NDEBUG */
alpar@9	2878
alpar@9	2879 Colamd_Col Col [], /* of size n_col+1 */
alpar@9	2880 Int A [], /* row indices of A */
alpar@9	2881 Int head [], /* head of degree lists and hash buckets */
alpar@9	2882 Int row_start, /* pointer to set of columns to check */
alpar@9	2883 Int row_length /* number of columns to check */
alpar@9	2884 )
alpar@9	2885 {
alpar@9	2886 /* === Local variables ================================================== */
alpar@9	2887
alpar@9	2888 Int hash ; /* hash value for a column */
alpar@9	2889 Int rp ; / pointer to a row */
alpar@9	2890 Int c ; /* a column index */
alpar@9	2891 Int super_c ; /* column index of the column to absorb into */
alpar@9	2892 Int cp1 ; / column pointer for column super_c */
alpar@9	2893 Int cp2 ; / column pointer for column c */
alpar@9	2894 Int length ; /* length of column super_c */
alpar@9	2895 Int prev_c ; /* column preceding c in hash bucket */
alpar@9	2896 Int i ; /* loop counter */
alpar@9	2897 Int rp_end ; / pointer to the end of the row */
alpar@9	2898 Int col ; /* a column index in the row to check */
alpar@9	2899 Int head_column ; /* first column in hash bucket or degree list */
alpar@9	2900 Int first_col ; /* first column in hash bucket */
alpar@9	2901
alpar@9	2902 /* === Consider each column in the row ================================== */
alpar@9	2903
alpar@9	2904 rp = &A [row_start] ;
alpar@9	2905 rp_end = rp + row_length ;
alpar@9	2906 while (rp < rp_end)
alpar@9	2907 {
alpar@9	2908 col = *rp++ ;
alpar@9	2909 if (COL_IS_DEAD (col))
alpar@9	2910 {
alpar@9	2911 continue ;
alpar@9	2912 }
alpar@9	2913
alpar@9	2914 /* get hash number for this column */
alpar@9	2915 hash = Col [col].shared3.hash ;
alpar@9	2916 ASSERT (hash <= n_col) ;
alpar@9	2917
alpar@9	2918 /* === Get the first column in this hash bucket ===================== */
alpar@9	2919
alpar@9	2920 head_column = head [hash] ;
alpar@9	2921 if (head_column > EMPTY)
alpar@9	2922 {
alpar@9	2923 first_col = Col [head_column].shared3.headhash ;
alpar@9	2924 }
alpar@9	2925 else
alpar@9	2926 {
alpar@9	2927 first_col = - (head_column + 2) ;
alpar@9	2928 }
alpar@9	2929
alpar@9	2930 /* === Consider each column in the hash bucket ====================== */
alpar@9	2931
alpar@9	2932 for (super_c = first_col ; super_c != EMPTY ;
alpar@9	2933 super_c = Col [super_c].shared4.hash_next)
alpar@9	2934 {
alpar@9	2935 ASSERT (COL_IS_ALIVE (super_c)) ;
alpar@9	2936 ASSERT (Col [super_c].shared3.hash == hash) ;
alpar@9	2937 length = Col [super_c].length ;
alpar@9	2938
alpar@9	2939 /* prev_c is the column preceding column c in the hash bucket */
alpar@9	2940 prev_c = super_c ;
alpar@9	2941
alpar@9	2942 /* === Compare super_c with all columns after it ================ */
alpar@9	2943
alpar@9	2944 for (c = Col [super_c].shared4.hash_next ;
alpar@9	2945 c != EMPTY ; c = Col [c].shared4.hash_next)
alpar@9	2946 {
alpar@9	2947 ASSERT (c != super_c) ;
alpar@9	2948 ASSERT (COL_IS_ALIVE (c)) ;
alpar@9	2949 ASSERT (Col [c].shared3.hash == hash) ;
alpar@9	2950
alpar@9	2951 /* not identical if lengths or scores are different */
alpar@9	2952 if (Col [c].length != length \|\|
alpar@9	2953 Col [c].shared2.score != Col [super_c].shared2.score)
alpar@9	2954 {
alpar@9	2955 prev_c = c ;
alpar@9	2956 continue ;
alpar@9	2957 }
alpar@9	2958
alpar@9	2959 /* compare the two columns */
alpar@9	2960 cp1 = &A [Col [super_c].start] ;
alpar@9	2961 cp2 = &A [Col [c].start] ;
alpar@9	2962
alpar@9	2963 for (i = 0 ; i < length ; i++)
alpar@9	2964 {
alpar@9	2965 /* the columns are "clean" (no dead rows) */
alpar@9	2966 ASSERT (ROW_IS_ALIVE (*cp1)) ;
alpar@9	2967 ASSERT (ROW_IS_ALIVE (*cp2)) ;
alpar@9	2968 /* row indices will same order for both supercols, */
alpar@9	2969 /* no gather scatter nessasary */
alpar@9	2970 if (cp1++ != cp2++)
alpar@9	2971 {
alpar@9	2972 break ;
alpar@9	2973 }
alpar@9	2974 }
alpar@9	2975
alpar@9	2976 /* the two columns are different if the for-loop "broke" */
alpar@9	2977 if (i != length)
alpar@9	2978 {
alpar@9	2979 prev_c = c ;
alpar@9	2980 continue ;
alpar@9	2981 }
alpar@9	2982
alpar@9	2983 /* === Got it! two columns are identical =================== */
alpar@9	2984
alpar@9	2985 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
alpar@9	2986
alpar@9	2987 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
alpar@9	2988 Col [c].shared1.parent = super_c ;
alpar@9	2989 KILL_NON_PRINCIPAL_COL (c) ;
alpar@9	2990 /* order c later, in order_children() */
alpar@9	2991 Col [c].shared2.order = EMPTY ;
alpar@9	2992 /* remove c from hash bucket */
alpar@9	2993 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
alpar@9	2994 }
alpar@9	2995 }
alpar@9	2996
alpar@9	2997 /* === Empty this hash bucket ======================================= */
alpar@9	2998
alpar@9	2999 if (head_column > EMPTY)
alpar@9	3000 {
alpar@9	3001 /* corresponding degree list "hash" is not empty */
alpar@9	3002 Col [head_column].shared3.headhash = EMPTY ;
alpar@9	3003 }
alpar@9	3004 else
alpar@9	3005 {
alpar@9	3006 /* corresponding degree list "hash" is empty */
alpar@9	3007 head [hash] = EMPTY ;
alpar@9	3008 }
alpar@9	3009 }
alpar@9	3010 }
alpar@9	3011
alpar@9	3012
alpar@9	3013 /* ========================================================================== */
alpar@9	3014 /* === garbage_collection =================================================== */
alpar@9	3015 /* ========================================================================== */
alpar@9	3016
alpar@9	3017 /*
alpar@9	3018 Defragments and compacts columns and rows in the workspace A. Used when
alpar@9	3019 all avaliable memory has been used while performing row merging. Returns
alpar@9	3020 the index of the first free position in A, after garbage collection. The
alpar@9	3021 time taken by this routine is linear is the size of the array A, which is
alpar@9	3022 itself linear in the number of nonzeros in the input matrix.
alpar@9	3023 Not user-callable.
alpar@9	3024 */
alpar@9	3025
alpar@9	3026 PRIVATE Int garbage_collection /* returns the new value of pfree */
alpar@9	3027 (
alpar@9	3028 /* === Parameters ======================================================= */
alpar@9	3029
alpar@9	3030 Int n_row, /* number of rows */
alpar@9	3031 Int n_col, /* number of columns */
alpar@9	3032 Colamd_Row Row [], /* row info */
alpar@9	3033 Colamd_Col Col [], /* column info */
alpar@9	3034 Int A [], /* A [0 ... Alen-1] holds the matrix */
alpar@9	3035 Int pfree / &A [0] ... pfree is in use */
alpar@9	3036 )
alpar@9	3037 {
alpar@9	3038 /* === Local variables ================================================== */
alpar@9	3039
alpar@9	3040 Int psrc ; / source pointer */
alpar@9	3041 Int pdest ; / destination pointer */
alpar@9	3042 Int j ; /* counter */
alpar@9	3043 Int r ; /* a row index */
alpar@9	3044 Int c ; /* a column index */
alpar@9	3045 Int length ; /* length of a row or column */
alpar@9	3046
alpar@9	3047 #ifndef NDEBUG
alpar@9	3048 Int debug_rows ;
alpar@9	3049 DEBUG2 (("Defrag..\n")) ;
alpar@9	3050 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
alpar@9	3051 debug_rows = 0 ;
alpar@9	3052 #endif /* NDEBUG */
alpar@9	3053
alpar@9	3054 /* === Defragment the columns =========================================== */
alpar@9	3055
alpar@9	3056 pdest = &A[0] ;
alpar@9	3057 for (c = 0 ; c < n_col ; c++)
alpar@9	3058 {
alpar@9	3059 if (COL_IS_ALIVE (c))
alpar@9	3060 {
alpar@9	3061 psrc = &A [Col [c].start] ;
alpar@9	3062
alpar@9	3063 /* move and compact the column */
alpar@9	3064 ASSERT (pdest <= psrc) ;
alpar@9	3065 Col [c].start = (Int) (pdest - &A [0]) ;
alpar@9	3066 length = Col [c].length ;
alpar@9	3067 for (j = 0 ; j < length ; j++)
alpar@9	3068 {
alpar@9	3069 r = *psrc++ ;
alpar@9	3070 if (ROW_IS_ALIVE (r))
alpar@9	3071 {
alpar@9	3072 *pdest++ = r ;
alpar@9	3073 }
alpar@9	3074 }
alpar@9	3075 Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
alpar@9	3076 }
alpar@9	3077 }
alpar@9	3078
alpar@9	3079 /* === Prepare to defragment the rows =================================== */
alpar@9	3080
alpar@9	3081 for (r = 0 ; r < n_row ; r++)
alpar@9	3082 {
alpar@9	3083 if (ROW_IS_DEAD (r) \|\| (Row [r].length == 0))
alpar@9	3084 {
alpar@9	3085 /* This row is already dead, or is of zero length. Cannot compact
alpar@9	3086 * a row of zero length, so kill it. NOTE: in the current version,
alpar@9	3087 * there are no zero-length live rows. Kill the row (for the first
alpar@9	3088 * time, or again) just to be safe. */
alpar@9	3089 KILL_ROW (r) ;
alpar@9	3090 }
alpar@9	3091 else
alpar@9	3092 {
alpar@9	3093 /* save first column index in Row [r].shared2.first_column */
alpar@9	3094 psrc = &A [Row [r].start] ;
alpar@9	3095 Row [r].shared2.first_column = *psrc ;
alpar@9	3096 ASSERT (ROW_IS_ALIVE (r)) ;
alpar@9	3097 /* flag the start of the row with the one's complement of row */
alpar@9	3098 *psrc = ONES_COMPLEMENT (r) ;
alpar@9	3099 #ifndef NDEBUG
alpar@9	3100 debug_rows++ ;
alpar@9	3101 #endif /* NDEBUG */
alpar@9	3102 }
alpar@9	3103 }
alpar@9	3104
alpar@9	3105 /* === Defragment the rows ============================================== */
alpar@9	3106
alpar@9	3107 psrc = pdest ;
alpar@9	3108 while (psrc < pfree)
alpar@9	3109 {
alpar@9	3110 /* find a negative number ... the start of a row */
alpar@9	3111 if (*psrc++ < 0)
alpar@9	3112 {
alpar@9	3113 psrc-- ;
alpar@9	3114 /* get the row index */
alpar@9	3115 r = ONES_COMPLEMENT (*psrc) ;
alpar@9	3116 ASSERT (r >= 0 && r < n_row) ;
alpar@9	3117 /* restore first column index */
alpar@9	3118 *psrc = Row [r].shared2.first_column ;
alpar@9	3119 ASSERT (ROW_IS_ALIVE (r)) ;
alpar@9	3120 ASSERT (Row [r].length > 0) ;
alpar@9	3121 /* move and compact the row */
alpar@9	3122 ASSERT (pdest <= psrc) ;
alpar@9	3123 Row [r].start = (Int) (pdest - &A [0]) ;
alpar@9	3124 length = Row [r].length ;
alpar@9	3125 for (j = 0 ; j < length ; j++)
alpar@9	3126 {
alpar@9	3127 c = *psrc++ ;
alpar@9	3128 if (COL_IS_ALIVE (c))
alpar@9	3129 {
alpar@9	3130 *pdest++ = c ;
alpar@9	3131 }
alpar@9	3132 }
alpar@9	3133 Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
alpar@9	3134 ASSERT (Row [r].length > 0) ;
alpar@9	3135 #ifndef NDEBUG
alpar@9	3136 debug_rows-- ;
alpar@9	3137 #endif /* NDEBUG */
alpar@9	3138 }
alpar@9	3139 }
alpar@9	3140 /* ensure we found all the rows */
alpar@9	3141 ASSERT (debug_rows == 0) ;
alpar@9	3142
alpar@9	3143 /* === Return the new value of pfree ==================================== */
alpar@9	3144
alpar@9	3145 return ((Int) (pdest - &A [0])) ;
alpar@9	3146 }
alpar@9	3147
alpar@9	3148
alpar@9	3149 /* ========================================================================== */
alpar@9	3150 /* === clear_mark =========================================================== */
alpar@9	3151 /* ========================================================================== */
alpar@9	3152
alpar@9	3153 /*
alpar@9	3154 Clears the Row [].shared2.mark array, and returns the new tag_mark.
alpar@9	3155 Return value is the new tag_mark. Not user-callable.
alpar@9	3156 */
alpar@9	3157
alpar@9	3158 PRIVATE Int clear_mark /* return the new value for tag_mark */
alpar@9	3159 (
alpar@9	3160 /* === Parameters ======================================================= */
alpar@9	3161
alpar@9	3162 Int tag_mark, /* new value of tag_mark */
alpar@9	3163 Int max_mark, /* max allowed value of tag_mark */
alpar@9	3164
alpar@9	3165 Int n_row, /* number of rows in A */
alpar@9	3166 Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
alpar@9	3167 )
alpar@9	3168 {
alpar@9	3169 /* === Local variables ================================================== */
alpar@9	3170
alpar@9	3171 Int r ;
alpar@9	3172
alpar@9	3173 if (tag_mark <= 0 \|\| tag_mark >= max_mark)
alpar@9	3174 {
alpar@9	3175 for (r = 0 ; r < n_row ; r++)
alpar@9	3176 {
alpar@9	3177 if (ROW_IS_ALIVE (r))
alpar@9	3178 {
alpar@9	3179 Row [r].shared2.mark = 0 ;
alpar@9	3180 }
alpar@9	3181 }
alpar@9	3182 tag_mark = 1 ;
alpar@9	3183 }
alpar@9	3184
alpar@9	3185 return (tag_mark) ;
alpar@9	3186 }
alpar@9	3187
alpar@9	3188
alpar@9	3189 /* ========================================================================== */
alpar@9	3190 /* === print_report ========================================================= */
alpar@9	3191 /* ========================================================================== */
alpar@9	3192
alpar@9	3193 PRIVATE void print_report
alpar@9	3194 (
alpar@9	3195 char *method,
alpar@9	3196 Int stats [COLAMD_STATS]
alpar@9	3197 )
alpar@9	3198 {
alpar@9	3199
alpar@9	3200 Int i1, i2, i3 ;
alpar@9	3201
alpar@9	3202 PRINTF (("\n%s version %d.%d, %s: ", method,
alpar@9	3203 COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
alpar@9	3204
alpar@9	3205 if (!stats)
alpar@9	3206 {
alpar@9	3207 PRINTF (("No statistics available.\n")) ;
alpar@9	3208 return ;
alpar@9	3209 }
alpar@9	3210
alpar@9	3211 i1 = stats [COLAMD_INFO1] ;
alpar@9	3212 i2 = stats [COLAMD_INFO2] ;
alpar@9	3213 i3 = stats [COLAMD_INFO3] ;
alpar@9	3214
alpar@9	3215 if (stats [COLAMD_STATUS] >= 0)
alpar@9	3216 {
alpar@9	3217 PRINTF (("OK. ")) ;
alpar@9	3218 }
alpar@9	3219 else
alpar@9	3220 {
alpar@9	3221 PRINTF (("ERROR. ")) ;
alpar@9	3222 }
alpar@9	3223
alpar@9	3224 switch (stats [COLAMD_STATUS])
alpar@9	3225 {
alpar@9	3226
alpar@9	3227 case COLAMD_OK_BUT_JUMBLED:
alpar@9	3228
alpar@9	3229 PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
alpar@9	3230
alpar@9	3231 PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
alpar@9	3232 method, i3)) ;
alpar@9	3233
alpar@9	3234 PRINTF(("%s: last seen duplicate or out-of-order row index: %d\n",
alpar@9	3235 method, INDEX (i2))) ;
alpar@9	3236
alpar@9	3237 PRINTF(("%s: last seen in column: %d",
alpar@9	3238 method, INDEX (i1))) ;
alpar@9	3239
alpar@9	3240 /* no break - fall through to next case instead */
alpar@9	3241
alpar@9	3242 case COLAMD_OK:
alpar@9	3243
alpar@9	3244 PRINTF(("\n")) ;
alpar@9	3245
alpar@9	3246 PRINTF(("%s: number of dense or empty rows ignored: %d\n",
alpar@9	3247 method, stats [COLAMD_DENSE_ROW])) ;
alpar@9	3248
alpar@9	3249 PRINTF(("%s: number of dense or empty columns ignored: %d\n",
alpar@9	3250 method, stats [COLAMD_DENSE_COL])) ;
alpar@9	3251
alpar@9	3252 PRINTF(("%s: number of garbage collections performed: %d\n",
alpar@9	3253 method, stats [COLAMD_DEFRAG_COUNT])) ;
alpar@9	3254 break ;
alpar@9	3255
alpar@9	3256 case COLAMD_ERROR_A_not_present:
alpar@9	3257
alpar@9	3258 PRINTF(("Array A (row indices of matrix) not present.\n")) ;
alpar@9	3259 break ;
alpar@9	3260
alpar@9	3261 case COLAMD_ERROR_p_not_present:
alpar@9	3262
alpar@9	3263 PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
alpar@9	3264 break ;
alpar@9	3265
alpar@9	3266 case COLAMD_ERROR_nrow_negative:
alpar@9	3267
alpar@9	3268 PRINTF(("Invalid number of rows (%d).\n", i1)) ;
alpar@9	3269 break ;
alpar@9	3270
alpar@9	3271 case COLAMD_ERROR_ncol_negative:
alpar@9	3272
alpar@9	3273 PRINTF(("Invalid number of columns (%d).\n", i1)) ;
alpar@9	3274 break ;
alpar@9	3275
alpar@9	3276 case COLAMD_ERROR_nnz_negative:
alpar@9	3277
alpar@9	3278 PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
alpar@9	3279 break ;
alpar@9	3280
alpar@9	3281 case COLAMD_ERROR_p0_nonzero:
alpar@9	3282
alpar@9	3283 PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
alpar@9	3284 break ;
alpar@9	3285
alpar@9	3286 case COLAMD_ERROR_A_too_small:
alpar@9	3287
alpar@9	3288 PRINTF(("Array A too small.\n")) ;
alpar@9	3289 PRINTF((" Need Alen >= %d, but given only Alen = %d.\n",
alpar@9	3290 i1, i2)) ;
alpar@9	3291 break ;
alpar@9	3292
alpar@9	3293 case COLAMD_ERROR_col_length_negative:
alpar@9	3294
alpar@9	3295 PRINTF
alpar@9	3296 (("Column %d has a negative number of nonzero entries (%d).\n",
alpar@9	3297 INDEX (i1), i2)) ;
alpar@9	3298 break ;
alpar@9	3299
alpar@9	3300 case COLAMD_ERROR_row_index_out_of_bounds:
alpar@9	3301
alpar@9	3302 PRINTF
alpar@9	3303 (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
alpar@9	3304 INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
alpar@9	3305 break ;
alpar@9	3306
alpar@9	3307 case COLAMD_ERROR_out_of_memory:
alpar@9	3308
alpar@9	3309 PRINTF(("Out of memory.\n")) ;
alpar@9	3310 break ;
alpar@9	3311
alpar@9	3312 /* v2.4: internal-error case deleted */
alpar@9	3313 }
alpar@9	3314 }
alpar@9	3315
alpar@9	3316
alpar@9	3317
alpar@9	3318
alpar@9	3319 /* ========================================================================== */
alpar@9	3320 /* === colamd debugging routines ============================================ */
alpar@9	3321 /* ========================================================================== */
alpar@9	3322
alpar@9	3323 /* When debugging is disabled, the remainder of this file is ignored. */
alpar@9	3324
alpar@9	3325 #ifndef NDEBUG
alpar@9	3326
alpar@9	3327
alpar@9	3328 /* ========================================================================== */
alpar@9	3329 /* === debug_structures ===================================================== */
alpar@9	3330 /* ========================================================================== */
alpar@9	3331
alpar@9	3332 /*
alpar@9	3333 At this point, all empty rows and columns are dead. All live columns
alpar@9	3334 are "clean" (containing no dead rows) and simplicial (no supercolumns
alpar@9	3335 yet). Rows may contain dead columns, but all live rows contain at
alpar@9	3336 least one live column.
alpar@9	3337 */
alpar@9	3338
alpar@9	3339 PRIVATE void debug_structures
alpar@9	3340 (
alpar@9	3341 /* === Parameters ======================================================= */
alpar@9	3342
alpar@9	3343 Int n_row,
alpar@9	3344 Int n_col,
alpar@9	3345 Colamd_Row Row [],
alpar@9	3346 Colamd_Col Col [],
alpar@9	3347 Int A [],
alpar@9	3348 Int n_col2
alpar@9	3349 )
alpar@9	3350 {
alpar@9	3351 /* === Local variables ================================================== */
alpar@9	3352
alpar@9	3353 Int i ;
alpar@9	3354 Int c ;
alpar@9	3355 Int *cp ;
alpar@9	3356 Int *cp_end ;
alpar@9	3357 Int len ;
alpar@9	3358 Int score ;
alpar@9	3359 Int r ;
alpar@9	3360 Int *rp ;
alpar@9	3361 Int *rp_end ;
alpar@9	3362 Int deg ;
alpar@9	3363
alpar@9	3364 /* === Check A, Row, and Col ============================================ */
alpar@9	3365
alpar@9	3366 for (c = 0 ; c < n_col ; c++)
alpar@9	3367 {
alpar@9	3368 if (COL_IS_ALIVE (c))
alpar@9	3369 {
alpar@9	3370 len = Col [c].length ;
alpar@9	3371 score = Col [c].shared2.score ;
alpar@9	3372 DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
alpar@9	3373 ASSERT (len > 0) ;
alpar@9	3374 ASSERT (score >= 0) ;
alpar@9	3375 ASSERT (Col [c].shared1.thickness == 1) ;
alpar@9	3376 cp = &A [Col [c].start] ;
alpar@9	3377 cp_end = cp + len ;
alpar@9	3378 while (cp < cp_end)
alpar@9	3379 {
alpar@9	3380 r = *cp++ ;
alpar@9	3381 ASSERT (ROW_IS_ALIVE (r)) ;
alpar@9	3382 }
alpar@9	3383 }
alpar@9	3384 else
alpar@9	3385 {
alpar@9	3386 i = Col [c].shared2.order ;
alpar@9	3387 ASSERT (i >= n_col2 && i < n_col) ;
alpar@9	3388 }
alpar@9	3389 }
alpar@9	3390
alpar@9	3391 for (r = 0 ; r < n_row ; r++)
alpar@9	3392 {
alpar@9	3393 if (ROW_IS_ALIVE (r))
alpar@9	3394 {
alpar@9	3395 i = 0 ;
alpar@9	3396 len = Row [r].length ;
alpar@9	3397 deg = Row [r].shared1.degree ;
alpar@9	3398 ASSERT (len > 0) ;
alpar@9	3399 ASSERT (deg > 0) ;
alpar@9	3400 rp = &A [Row [r].start] ;
alpar@9	3401 rp_end = rp + len ;
alpar@9	3402 while (rp < rp_end)
alpar@9	3403 {
alpar@9	3404 c = *rp++ ;
alpar@9	3405 if (COL_IS_ALIVE (c))
alpar@9	3406 {
alpar@9	3407 i++ ;
alpar@9	3408 }
alpar@9	3409 }
alpar@9	3410 ASSERT (i > 0) ;
alpar@9	3411 }
alpar@9	3412 }
alpar@9	3413 }
alpar@9	3414
alpar@9	3415
alpar@9	3416 /* ========================================================================== */
alpar@9	3417 /* === debug_deg_lists ====================================================== */
alpar@9	3418 /* ========================================================================== */
alpar@9	3419
alpar@9	3420 /*
alpar@9	3421 Prints the contents of the degree lists. Counts the number of columns
alpar@9	3422 in the degree list and compares it to the total it should have. Also
alpar@9	3423 checks the row degrees.
alpar@9	3424 */
alpar@9	3425
alpar@9	3426 PRIVATE void debug_deg_lists
alpar@9	3427 (
alpar@9	3428 /* === Parameters ======================================================= */
alpar@9	3429
alpar@9	3430 Int n_row,
alpar@9	3431 Int n_col,
alpar@9	3432 Colamd_Row Row [],
alpar@9	3433 Colamd_Col Col [],
alpar@9	3434 Int head [],
alpar@9	3435 Int min_score,
alpar@9	3436 Int should,
alpar@9	3437 Int max_deg
alpar@9	3438 )
alpar@9	3439 {
alpar@9	3440 /* === Local variables ================================================== */
alpar@9	3441
alpar@9	3442 Int deg ;
alpar@9	3443 Int col ;
alpar@9	3444 Int have ;
alpar@9	3445 Int row ;
alpar@9	3446
alpar@9	3447 /* === Check the degree lists =========================================== */
alpar@9	3448
alpar@9	3449 if (n_col > 10000 && colamd_debug <= 0)
alpar@9	3450 {
alpar@9	3451 return ;
alpar@9	3452 }
alpar@9	3453 have = 0 ;
alpar@9	3454 DEBUG4 (("Degree lists: %d\n", min_score)) ;
alpar@9	3455 for (deg = 0 ; deg <= n_col ; deg++)
alpar@9	3456 {
alpar@9	3457 col = head [deg] ;
alpar@9	3458 if (col == EMPTY)
alpar@9	3459 {
alpar@9	3460 continue ;
alpar@9	3461 }
alpar@9	3462 DEBUG4 (("%d:", deg)) ;
alpar@9	3463 while (col != EMPTY)
alpar@9	3464 {
alpar@9	3465 DEBUG4 ((" %d", col)) ;
alpar@9	3466 have += Col [col].shared1.thickness ;
alpar@9	3467 ASSERT (COL_IS_ALIVE (col)) ;
alpar@9	3468 col = Col [col].shared4.degree_next ;
alpar@9	3469 }
alpar@9	3470 DEBUG4 (("\n")) ;
alpar@9	3471 }
alpar@9	3472 DEBUG4 (("should %d have %d\n", should, have)) ;
alpar@9	3473 ASSERT (should == have) ;
alpar@9	3474
alpar@9	3475 /* === Check the row degrees ============================================ */
alpar@9	3476
alpar@9	3477 if (n_row > 10000 && colamd_debug <= 0)
alpar@9	3478 {
alpar@9	3479 return ;
alpar@9	3480 }
alpar@9	3481 for (row = 0 ; row < n_row ; row++)
alpar@9	3482 {
alpar@9	3483 if (ROW_IS_ALIVE (row))
alpar@9	3484 {
alpar@9	3485 ASSERT (Row [row].shared1.degree <= max_deg) ;
alpar@9	3486 }
alpar@9	3487 }
alpar@9	3488 }
alpar@9	3489
alpar@9	3490
alpar@9	3491 /* ========================================================================== */
alpar@9	3492 /* === debug_mark =========================================================== */
alpar@9	3493 /* ========================================================================== */
alpar@9	3494
alpar@9	3495 /*
alpar@9	3496 Ensures that the tag_mark is less that the maximum and also ensures that
alpar@9	3497 each entry in the mark array is less than the tag mark.
alpar@9	3498 */
alpar@9	3499
alpar@9	3500 PRIVATE void debug_mark
alpar@9	3501 (
alpar@9	3502 /* === Parameters ======================================================= */
alpar@9	3503
alpar@9	3504 Int n_row,
alpar@9	3505 Colamd_Row Row [],
alpar@9	3506 Int tag_mark,
alpar@9	3507 Int max_mark
alpar@9	3508 )
alpar@9	3509 {
alpar@9	3510 /* === Local variables ================================================== */
alpar@9	3511
alpar@9	3512 Int r ;
alpar@9	3513
alpar@9	3514 /* === Check the Row marks ============================================== */
alpar@9	3515
alpar@9	3516 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
alpar@9	3517 if (n_row > 10000 && colamd_debug <= 0)
alpar@9	3518 {
alpar@9	3519 return ;
alpar@9	3520 }
alpar@9	3521 for (r = 0 ; r < n_row ; r++)
alpar@9	3522 {
alpar@9	3523 ASSERT (Row [r].shared2.mark < tag_mark) ;
alpar@9	3524 }
alpar@9	3525 }
alpar@9	3526
alpar@9	3527
alpar@9	3528 /* ========================================================================== */
alpar@9	3529 /* === debug_matrix ========================================================= */
alpar@9	3530 /* ========================================================================== */
alpar@9	3531
alpar@9	3532 /*
alpar@9	3533 Prints out the contents of the columns and the rows.
alpar@9	3534 */
alpar@9	3535
alpar@9	3536 PRIVATE void debug_matrix
alpar@9	3537 (
alpar@9	3538 /* === Parameters ======================================================= */
alpar@9	3539
alpar@9	3540 Int n_row,
alpar@9	3541 Int n_col,
alpar@9	3542 Colamd_Row Row [],
alpar@9	3543 Colamd_Col Col [],
alpar@9	3544 Int A []
alpar@9	3545 )
alpar@9	3546 {
alpar@9	3547 /* === Local variables ================================================== */
alpar@9	3548
alpar@9	3549 Int r ;
alpar@9	3550 Int c ;
alpar@9	3551 Int *rp ;
alpar@9	3552 Int *rp_end ;
alpar@9	3553 Int *cp ;
alpar@9	3554 Int *cp_end ;
alpar@9	3555
alpar@9	3556 /* === Dump the rows and columns of the matrix ========================== */
alpar@9	3557
alpar@9	3558 if (colamd_debug < 3)
alpar@9	3559 {
alpar@9	3560 return ;
alpar@9	3561 }
alpar@9	3562 DEBUG3 (("DUMP MATRIX:\n")) ;
alpar@9	3563 for (r = 0 ; r < n_row ; r++)
alpar@9	3564 {
alpar@9	3565 DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
alpar@9	3566 if (ROW_IS_DEAD (r))
alpar@9	3567 {
alpar@9	3568 continue ;
alpar@9	3569 }
alpar@9	3570 DEBUG3 (("start %d length %d degree %d\n",
alpar@9	3571 Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
alpar@9	3572 rp = &A [Row [r].start] ;
alpar@9	3573 rp_end = rp + Row [r].length ;
alpar@9	3574 while (rp < rp_end)
alpar@9	3575 {
alpar@9	3576 c = *rp++ ;
alpar@9	3577 DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ;
alpar@9	3578 }
alpar@9	3579 }
alpar@9	3580
alpar@9	3581 for (c = 0 ; c < n_col ; c++)
alpar@9	3582 {
alpar@9	3583 DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
alpar@9	3584 if (COL_IS_DEAD (c))
alpar@9	3585 {
alpar@9	3586 continue ;
alpar@9	3587 }
alpar@9	3588 DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
alpar@9	3589 Col [c].start, Col [c].length,
alpar@9	3590 Col [c].shared1.thickness, Col [c].shared2.score)) ;
alpar@9	3591 cp = &A [Col [c].start] ;
alpar@9	3592 cp_end = cp + Col [c].length ;
alpar@9	3593 while (cp < cp_end)
alpar@9	3594 {
alpar@9	3595 r = *cp++ ;
alpar@9	3596 DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ;
alpar@9	3597 }
alpar@9	3598 }
alpar@9	3599 }
alpar@9	3600
alpar@9	3601 PRIVATE void colamd_get_debug
alpar@9	3602 (
alpar@9	3603 char *method
alpar@9	3604 )
alpar@9	3605 {
alpar@9	3606 FILE *f ;
alpar@9	3607 colamd_debug = 0 ; /* no debug printing */
alpar@9	3608 f = fopen ("debug", "r") ;
alpar@9	3609 if (f == (FILE *) NULL)
alpar@9	3610 {
alpar@9	3611 colamd_debug = 0 ;
alpar@9	3612 }
alpar@9	3613 else
alpar@9	3614 {
alpar@9	3615 fscanf (f, "%d", &colamd_debug) ;
alpar@9	3616 fclose (f) ;
alpar@9	3617 }
alpar@9	3618 DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
alpar@9	3619 method, colamd_debug)) ;
alpar@9	3620 }
alpar@9	3621
alpar@9	3622 #endif /* NDEBUG */