Context Navigation

colamd.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f 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

Line
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 @ 1:c445c931472f

Download in other formats: