Context Navigation

colamd.c @ 2:4c8956a7bdf4

Last change on this file since 2:4c8956a7bdf4 was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

Generated files and doc/notes are removed

File size: 124.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: glpk-cmake/src/colamd/colamd.c @ 2:4c8956a7bdf4

Download in other formats: