Commits

Shlomi Fish committed c1cac07

Add more files.

Comments (0)

Files changed (14)

+/* ========================================================================== */
+/* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
+/* ========================================================================== */
+
+/*
+    colamd:  an approximate minimum degree column ordering algorithm,
+    	for LU factorization of symmetric or unsymmetric matrices,
+	QR factorization, least squares, interior point methods for
+	linear programming problems, and other related problems.
+
+    symamd:  an approximate minimum degree ordering algorithm for Cholesky
+    	factorization of symmetric matrices.
+
+    Purpose:
+
+	Colamd computes a permutation Q such that the Cholesky factorization of
+	(AQ)'(AQ) has less fill-in and requires fewer floating point operations
+	than A'A.  This also provides a good ordering for sparse partial
+	pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
+	factorization, and P is computed during numerical factorization via
+	conventional partial pivoting with row interchanges.  Colamd is the
+	column ordering method used in SuperLU, part of the ScaLAPACK library.
+	It is also available as built-in function in Matlab Version 6,
+	available from MathWorks, Inc. (http://www.mathworks.com).  This
+	routine can be used in place of colmmd in Matlab.  By default, the \
+	and / operators in Matlab perform a column ordering (using colmmd
+	or colamd) prior to LU factorization using sparse partial pivoting,
+	in the built-in Matlab lu(A) routine.
+
+    	Symamd computes a permutation P of a symmetric matrix A such that the
+	Cholesky factorization of PAP' has less fill-in and requires fewer
+	floating point operations than A.  Symamd constructs a matrix M such
+	that M'M has the same nonzero pattern of A, and then orders the columns
+	of M using colmmd.  The column ordering of M is then returned as the
+	row and column ordering P of A. 
+
+    Authors:
+
+	The authors of the code itself are Stefan I. Larimore and Timothy A.
+	Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was
+	developed in collaboration with John Gilbert, Xerox PARC, and Esmond
+	Ng, Oak Ridge National Laboratory.
+
+    Date:
+
+	May 4, 2001.  Version 2.1.
+
+    Acknowledgements:
+
+	This work was supported by the National Science Foundation, under
+	grants DMS-9504974 and DMS-9803599.
+
+    Notice:
+
+	Copyright (c) 1998-2001 by the University of Florida.
+	All Rights Reserved.
+
+	THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+	EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+
+	Permission is hereby granted to use or copy this program for any
+	purpose, provided the above notices are retained on all copies.
+	User documentation of any code that uses this code must cite the
+	Authors, the Copyright, and "Used by permission."  If this code is
+	accessible from within Matlab, then typing "help colamd" and "help
+	symamd" must cite the Authors.  Permission to modify the code and to
+	distribute modified code is granted, provided the above notices are
+	retained, and a notice that the code was modified is included with the
+	above copyright notice.  You must also retain the Availability
+	information below, of the original version.
+
+	This software is provided free of charge.
+
+    Availability:
+
+	The colamd/symamd library is available at
+
+	    http://www.cise.ufl.edu/research/sparse/colamd/
+
+	This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
+	file.  It requires the colamd.h file.  It is required by the colamdmex.c
+	and symamdmex.c files, for the Matlab interface to colamd and symamd.
+
+    Changes to the colamd library since Version 1.0 and 1.1:
+
+	No bugs were found in version 1.1.  These changes merely add new
+	functionality.
+
+    	* added the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro.
+
+	* moved the output statistics, from A, to a separate output argument.
+		The arguments changed for the C-callable routines.
+
+	* added colamd_report and symamd_report.
+
+	* added a C-callable symamd routine.  Formerly, symamd was only
+		available as a mexFunction from Matlab.
+
+	* added error-checking to symamd.  Formerly, it assumed its input
+		was error-free.
+
+	* added the optional stats and knobs arguments to the symamd mexFunction
+
+	* deleted colamd_help.  A help message is still available from
+		"help colamd" and "help symamd" in Matlab.
+
+	* deleted colamdtree.m and symamdtree.m.  Now, colamd.m and symamd.m
+		also do the elimination tree post-ordering.  The Version 1.1
+		colamd and symamd mexFunctions, which do not do the post-
+		ordering, are now visible as colamdmex and symamdmex from
+		Matlab.  Essentialy, the post-ordering is now the default
+		behavior of colamd.m and symamd.m, to match the behavior of
+		colmmd and symmmd.  The post-ordering is only available in the
+		Matlab interface, not the C-callable interface.
+
+	* made a slight change to the dense row/column detection in symamd,
+		to match the stated specifications.
+
+    Changes from Version 2.0 to 2.1:
+
+	* TRUE and FALSE are predefined on some systems, so they are defined
+		here only if not already defined.
+	
+	* web site changed
+
+	* UNIX Makefile modified, to handle the case if "." is not in your path.
+
+*/
+
+/* ========================================================================== */
+/* === Description of user-callable routines ================================ */
+/* ========================================================================== */
+
+/*
+    ----------------------------------------------------------------------------
+    colamd_recommended:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    int colamd_recommended (int nnz, int n_row, int n_col) ;
+
+	    or as a C macro
+
+	    #include "colamd.h"
+	    Alen = COLAMD_RECOMMENDED (int nnz, int n_row, int n_col) ;
+
+	Purpose:
+
+	    Returns recommended value of Alen for use by colamd.  Returns -1
+	    if any input argument is negative.  The use of this routine
+	    or macro is optional.  Note that the macro uses its arguments
+	    more than once, so be careful for side effects, if you pass
+	    expressions as arguments to COLAMD_RECOMMENDED.  Not needed for
+	    symamd, which dynamically allocates its own memory.
+
+	Arguments (all input arguments):
+
+	    int nnz ;		Number of nonzeros in the matrix A.  This must
+				be the same value as p [n_col] in the call to
+				colamd - otherwise you will get a wrong value
+				of the recommended memory to use.
+
+	    int n_row ;		Number of rows in the matrix A.
+
+	    int n_col ;		Number of columns in the matrix A.
+
+    ----------------------------------------------------------------------------
+    colamd_set_defaults:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    colamd_set_defaults (int knobs [COLAMD_KNOBS]) ;
+
+	Purpose:
+
+	    Sets the default parameters.  The use of this routine is optional.
+
+	Arguments:
+
+	    double knobs [COLAMD_KNOBS] ;	Output only.
+
+		Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col)
+		entries are removed prior to ordering.  Columns with more than
+		(knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to
+		ordering, and placed last in the output column ordering. 
+
+		Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
+		Rows and columns with more than (knobs [COLAMD_DENSE_ROW] * n)
+		entries are removed prior to ordering, and placed last in the
+		output ordering.
+
+		COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
+		respectively, in colamd.h.  Default values of these two knobs
+		are both 0.5.  Currently, only knobs [0] and knobs [1] are
+		used, but future versions may use more knobs.  If so, they will
+		be properly set to their defaults by the future version of
+		colamd_set_defaults, so that the code that calls colamd will
+		not need to change, assuming that you either use
+		colamd_set_defaults, or pass a (double *) NULL pointer as the
+		knobs array to colamd or symamd.
+
+    ----------------------------------------------------------------------------
+    colamd:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    int colamd (int n_row, int n_col, int Alen, int *A, int *p,
+	    	double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
+
+	Purpose:
+
+	    Computes a column ordering (Q) of A such that P(AQ)=LU or
+	    (AQ)'AQ=LL' have less fill-in and require fewer floating point
+	    operations than factorizing the unpermuted matrix A or A'A,
+	    respectively.
+	    
+	Returns:
+
+	    TRUE (1) if successful, FALSE (0) otherwise.
+
+	Arguments:
+
+	    int n_row ;		Input argument.
+
+		Number of rows in the matrix A.
+		Restriction:  n_row >= 0.
+		Colamd returns FALSE if n_row is negative.
+
+	    int n_col ;		Input argument.
+
+		Number of columns in the matrix A.
+		Restriction:  n_col >= 0.
+		Colamd returns FALSE if n_col is negative.
+
+	    int Alen ;		Input argument.
+
+		Restriction (see note):
+		Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
+		Colamd returns FALSE if these conditions are not met.
+
+		Note:  this restriction makes an modest assumption regarding
+		the size of the two typedef's structures in colamd.h.
+		We do, however, guarantee that
+
+			Alen >= colamd_recommended (nnz, n_row, n_col)
+		
+		or equivalently as a C preprocessor macro: 
+
+			Alen >= COLAMD_RECOMMENDED (nnz, n_row, n_col)
+
+		will be sufficient.
+
+	    int A [Alen] ;	Input argument, undefined on output.
+
+		A is an integer array of size Alen.  Alen must be at least as
+		large as the bare minimum value given above, but this is very
+		low, and can result in excessive run time.  For best
+		performance, we recommend that Alen be greater than or equal to
+		colamd_recommended (nnz, n_row, n_col), which adds
+		nnz/5 to the bare minimum value given above.
+
+		On input, the row indices of the entries in column c of the
+		matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
+		in a given column c need not be in ascending order, and
+		duplicate row indices may be be present.  However, colamd will
+		work a little faster if both of these conditions are met
+		(Colamd puts the matrix into this format, if it finds that the
+		the conditions are not met).
+
+		The matrix is 0-based.  That is, rows are in the range 0 to
+		n_row-1, and columns are in the range 0 to n_col-1.  Colamd
+		returns FALSE if any row index is out of range.
+
+		The contents of A are modified during ordering, and are
+		undefined on output.
+
+	    int p [n_col+1] ;	Both input and output argument.
+
+		p is an integer array of size n_col+1.  On input, it holds the
+		"pointers" for the column form of the matrix A.  Column c of
+		the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
+		entry, p [0], must be zero, and p [c] <= p [c+1] must hold
+		for all c in the range 0 to n_col-1.  The value p [n_col] is
+		thus the total number of entries in the pattern of the matrix A.
+		Colamd returns FALSE if these conditions are not met.
+
+		On output, if colamd returns TRUE, the array p holds the column
+		permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
+		the first column index in the new ordering, and p [n_col-1] is
+		the last.  That is, p [k] = j means that column j of A is the
+		kth pivot column, in AQ, where k is in the range 0 to n_col-1
+		(p [0] = j means that column j of A is the first column in AQ).
+
+		If colamd returns FALSE, then no permutation is returned, and
+		p is undefined on output.
+
+	    double knobs [COLAMD_KNOBS] ;	Input argument.
+
+		See colamd_set_defaults for a description.
+
+	    int stats [COLAMD_STATS] ;		Output argument.
+
+		Statistics on the ordering, and error status.
+		See colamd.h for related definitions.
+		Colamd returns FALSE if stats is not present.
+
+		stats [0]:  number of dense or empty rows ignored.
+
+		stats [1]:  number of dense or empty columns ignored (and
+				ordered last in the output permutation p)
+				Note that a row can become "empty" if it
+				contains only "dense" and/or "empty" columns,
+				and similarly a column can become "empty" if it
+				only contains "dense" and/or "empty" rows.
+
+		stats [2]:  number of garbage collections performed.
+				This can be excessively high if Alen is close
+				to the minimum required value.
+
+		stats [3]:  status code.  < 0 is an error code.
+			    > 1 is a warning or notice.
+
+			0	OK.  Each column of the input matrix contained
+				row indices in increasing order, with no
+				duplicates.
+
+			1	OK, but columns of input matrix were jumbled
+				(unsorted columns or duplicate entries).  Colamd
+				had to do some extra work to sort the matrix
+				first and remove duplicate entries, but it
+				still was able to return a valid permutation
+				(return value of colamd was TRUE).
+
+					stats [4]: highest numbered column that
+						is unsorted or has duplicate
+						entries.
+					stats [5]: last seen duplicate or
+						unsorted row index.
+					stats [6]: number of duplicate or
+						unsorted row indices.
+
+			-1	A is a null pointer
+
+			-2	p is a null pointer
+
+			-3 	n_row is negative
+
+					stats [4]: n_row
+
+			-4	n_col is negative
+
+					stats [4]: n_col
+
+			-5	number of nonzeros in matrix is negative
+
+					stats [4]: number of nonzeros, p [n_col]
+
+			-6	p [0] is nonzero
+
+					stats [4]: p [0]
+
+			-7	A is too small
+
+					stats [4]: required size
+					stats [5]: actual size (Alen)
+
+			-8	a column has a negative number of entries
+
+					stats [4]: column with < 0 entries
+					stats [5]: number of entries in col
+
+			-9	a row index is out of bounds
+
+					stats [4]: column with bad row index
+					stats [5]: bad row index
+					stats [6]: n_row, # of rows of matrx
+
+			-10	(unused; see symamd.c)
+
+			-999	(unused; see symamd.c)
+
+		Future versions may return more statistics in the stats array.
+
+	Example:
+	
+	    See http://www.cise.ufl.edu/research/sparse/colamd/example.c
+	    for a complete example.
+
+	    To order the columns of a 5-by-4 matrix with 11 nonzero entries in
+	    the following nonzero pattern
+
+	    	x 0 x 0
+		x 0 x x
+		0 x x 0
+		0 0 x x
+		x x 0 0
+
+	    with default knobs and no output statistics, do the following:
+
+		#include "colamd.h"
+		#define ALEN COLAMD_RECOMMENDED (11, 5, 4)
+		int A [ALEN] = {1, 2, 5, 3, 5, 1, 2, 3, 4, 2, 4} ;
+		int p [ ] = {0, 3, 5, 9, 11} ;
+		int stats [COLAMD_STATS] ;
+		colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
+
+	    The permutation is returned in the array p, and A is destroyed.
+
+    ----------------------------------------------------------------------------
+    symamd:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    int symamd (int n, int *A, int *p, int *perm,
+	    	int knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
+		void (*allocate) (size_t, size_t), void (*release) (void *)) ;
+
+	Purpose:
+
+    	    The symamd routine computes an ordering P of a symmetric sparse
+	    matrix A such that the Cholesky factorization PAP' = LL' remains
+	    sparse.  It is based on a column ordering of a matrix M constructed
+	    so that the nonzero pattern of M'M is the same as A.  The matrix A
+	    is assumed to be symmetric; only the strictly lower triangular part
+	    is accessed.  You must pass your selected memory allocator (usually
+	    calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
+	    memory for the temporary matrix M.
+
+	Returns:
+
+	    TRUE (1) if successful, FALSE (0) otherwise.
+
+	Arguments:
+
+	    int n ;		Input argument.
+
+	    	Number of rows and columns in the symmetrix matrix A.
+		Restriction:  n >= 0.
+		Symamd returns FALSE if n is negative.
+
+	    int A [nnz] ;	Input argument.
+
+	    	A is an integer array of size nnz, where nnz = p [n].
+		
+		The row indices of the entries in column c of the matrix are
+		held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
+		given column c need not be in ascending order, and duplicate
+		row indices may be present.  However, symamd will run faster
+		if the columns are in sorted order with no duplicate entries. 
+
+		The matrix is 0-based.  That is, rows are in the range 0 to
+		n-1, and columns are in the range 0 to n-1.  Symamd
+		returns FALSE if any row index is out of range.
+
+		The contents of A are not modified.
+
+	    int p [n+1] ;   	Input argument.
+
+		p is an integer array of size n+1.  On input, it holds the
+		"pointers" for the column form of the matrix A.  Column c of
+		the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
+		entry, p [0], must be zero, and p [c] <= p [c+1] must hold
+		for all c in the range 0 to n-1.  The value p [n] is
+		thus the total number of entries in the pattern of the matrix A.
+		Symamd returns FALSE if these conditions are not met.
+
+		The contents of p are not modified.
+
+	    int perm [n+1] ;   	Output argument.
+
+		On output, if symamd returns TRUE, the array perm holds the
+		permutation P, where perm [0] is the first index in the new
+		ordering, and perm [n-1] is the last.  That is, perm [k] = j
+		means that row and column j of A is the kth column in PAP',
+		where k is in the range 0 to n-1 (perm [0] = j means
+		that row and column j of A are the first row and column in
+		PAP').  The array is used as a workspace during the ordering,
+		which is why it must be of length n+1, not just n.
+
+	    double knobs [COLAMD_KNOBS] ;	Input argument.
+
+		See colamd_set_defaults for a description.
+
+	    int stats [COLAMD_STATS] ;		Output argument.
+
+		Statistics on the ordering, and error status.
+		See colamd.h for related definitions.
+		Symamd returns FALSE if stats is not present.
+
+		stats [0]:  number of dense or empty row and columns ignored
+				(and ordered last in the output permutation 
+				perm).  Note that a row/column can become
+				"empty" if it contains only "dense" and/or
+				"empty" columns/rows.
+
+		stats [1]:  (same as stats [0])
+
+		stats [2]:  number of garbage collections performed.
+
+		stats [3]:  status code.  < 0 is an error code.
+			    > 1 is a warning or notice.
+
+			0	OK.  Each column of the input matrix contained
+				row indices in increasing order, with no
+				duplicates.
+
+			1	OK, but columns of input matrix were jumbled
+				(unsorted columns or duplicate entries).  Symamd
+				had to do some extra work to sort the matrix
+				first and remove duplicate entries, but it
+				still was able to return a valid permutation
+				(return value of symamd was TRUE).
+
+					stats [4]: highest numbered column that
+						is unsorted or has duplicate
+						entries.
+					stats [5]: last seen duplicate or
+						unsorted row index.
+					stats [6]: number of duplicate or
+						unsorted row indices.
+
+			-1	A is a null pointer
+
+			-2	p is a null pointer
+
+			-3	(unused, see colamd.c)
+
+			-4 	n is negative
+
+					stats [4]: n
+
+			-5	number of nonzeros in matrix is negative
+
+					stats [4]: # of nonzeros (p [n]).
+
+			-6	p [0] is nonzero
+
+					stats [4]: p [0]
+
+			-7	(unused)
+
+			-8	a column has a negative number of entries
+
+					stats [4]: column with < 0 entries
+					stats [5]: number of entries in col
+
+			-9	a row index is out of bounds
+
+					stats [4]: column with bad row index
+					stats [5]: bad row index
+					stats [6]: n_row, # of rows of matrx
+
+			-10	out of memory (unable to allocate temporary
+				workspace for M or count arrays using the
+				"allocate" routine passed into symamd).
+
+			-999	internal error.  colamd failed to order the
+				matrix M, when it should have succeeded.  This
+				indicates a bug.  If this (and *only* this)
+				error code occurs, please contact the authors.
+				Don't contact the authors if you get any other
+				error code.
+
+		Future versions may return more statistics in the stats array.
+
+	    void * (*allocate) (size_t, size_t)
+
+	    	A pointer to a function providing memory allocation.  The
+		allocated memory must be returned initialized to zero.  For a
+		C application, this argument should normally be a pointer to
+		calloc.  For a Matlab mexFunction, the routine mxCalloc is
+		passed instead.
+
+	    void (*release) (size_t, size_t)
+
+	    	A pointer to a function that frees memory allocated by the
+		memory allocation routine above.  For a C application, this
+		argument should normally be a pointer to free.  For a Matlab
+		mexFunction, the routine mxFree is passed instead.
+
+
+    ----------------------------------------------------------------------------
+    colamd_report:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    colamd_report (int stats [COLAMD_STATS]) ;
+
+	Purpose:
+
+	    Prints the error status and statistics recorded in the stats
+	    array on the standard error output (for a standard C routine)
+	    or on the Matlab output (for a mexFunction).
+
+	Arguments:
+
+	    int stats [COLAMD_STATS] ;	Input only.  Statistics from colamd.
+
+
+    ----------------------------------------------------------------------------
+    symamd_report:
+    ----------------------------------------------------------------------------
+
+	C syntax:
+
+	    #include "colamd.h"
+	    symamd_report (int stats [COLAMD_STATS]) ;
+
+	Purpose:
+
+	    Prints the error status and statistics recorded in the stats
+	    array on the standard error output (for a standard C routine)
+	    or on the Matlab output (for a mexFunction).
+
+	Arguments:
+
+	    int stats [COLAMD_STATS] ;	Input only.  Statistics from symamd.
+
+
+*/
+
+/* ========================================================================== */
+/* === Scaffolding code definitions  ======================================== */
+/* ========================================================================== */
+
+/* Ensure that debugging is turned off: */
+#ifndef NDEBUG
+#define NDEBUG
+#endif /* NDEBUG */
+
+/*
+   Our "scaffolding code" philosophy:  In our opinion, well-written library
+   code should keep its "debugging" code, and just normally have it turned off
+   by the compiler so as not to interfere with performance.  This serves
+   several purposes:
+
+   (1) assertions act as comments to the reader, telling you what the code
+	expects at that point.  All assertions will always be true (unless
+	there really is a bug, of course).
+
+   (2) leaving in the scaffolding code assists anyone who would like to modify
+	the code, or understand the algorithm (by reading the debugging output,
+	one can get a glimpse into what the code is doing).
+
+   (3) (gasp!) for actually finding bugs.  This code has been heavily tested
+	and "should" be fully functional and bug-free ... but you never know...
+
+    To enable debugging, comment out the "#define NDEBUG" above.  For a Matlab
+    mexFunction, you will also need to modify mexopts.sh to remove the -DNDEBUG
+    definition.  The code will become outrageously slow when debugging is
+    enabled.  To control the level of debugging output, set an environment
+    variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
+    you should see the following message on the standard output:
+
+    	colamd: debug version, D = 1 (THIS WILL BE SLOW!)
+
+    or a similar message for symamd.  If you don't, then debugging has not
+    been enabled.
+
+*/
+
+/* ========================================================================== */
+/* === Include files ======================================================== */
+/* ========================================================================== */
+
+#include "colamd.h"
+#include <limits.h>
+
+#ifdef MATLAB_MEX_FILE
+#include "mex.h"
+#include "matrix.h"
+#else
+#include <stdio.h>
+#include <assert.h>
+#endif /* MATLAB_MEX_FILE */
+
+#ifdef FORTIFY
+# include "lp_fortify.h"
+#endif
+
+
+/* ========================================================================== */
+/* === Definitions ========================================================== */
+/* ========================================================================== */
+
+/* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
+#define PUBLIC
+#define PRIVATE static
+
+#define MAX(a,b) (((a) > (b)) ? (a) : (b))
+#define MIN(a,b) (((a) < (b)) ? (a) : (b))
+
+#define ONES_COMPLEMENT(r) (-(r)-1)
+
+/* -------------------------------------------------------------------------- */
+/* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
+/* -------------------------------------------------------------------------- */
+
+#ifndef TRUE
+#define TRUE (1)
+#endif
+
+#ifndef FALSE
+#define FALSE (0)
+#endif
+
+/* -------------------------------------------------------------------------- */
+
+#define EMPTY	(-1)
+
+/* Row and column status */
+#define ALIVE	(0)
+#define DEAD	(-1)
+
+/* Column status */
+#define DEAD_PRINCIPAL		(-1)
+#define DEAD_NON_PRINCIPAL	(-2)
+
+/* Macros for row and column status update and checking. */
+#define ROW_IS_DEAD(r)			ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
+#define ROW_IS_MARKED_DEAD(row_mark)	(row_mark < ALIVE)
+#define ROW_IS_ALIVE(r)			(Row [r].shared2.mark >= ALIVE)
+#define COL_IS_DEAD(c)			(Col [c].start < ALIVE)
+#define COL_IS_ALIVE(c)			(Col [c].start >= ALIVE)
+#define COL_IS_DEAD_PRINCIPAL(c)	(Col [c].start == DEAD_PRINCIPAL)
+#define KILL_ROW(r)			{ Row [r].shared2.mark = DEAD ; }
+#define KILL_PRINCIPAL_COL(c)		{ Col [c].start = DEAD_PRINCIPAL ; }
+#define KILL_NON_PRINCIPAL_COL(c)	{ Col [c].start = DEAD_NON_PRINCIPAL ; }
+
+/* ========================================================================== */
+/* === Colamd reporting mechanism =========================================== */
+/* ========================================================================== */
+
+#ifdef MATLAB_MEX_FILE
+
+/* use mexPrintf in a Matlab mexFunction, for debugging and statistics output */
+#define PRINTF mexPrintf
+
+/* In Matlab, matrices are 1-based to the user, but 0-based internally */
+#define INDEX(i) ((i)+1)
+
+#else
+
+/* Use printf in standard C environment, for debugging and statistics output. */
+/* Output is generated only if debugging is enabled at compile time, or if */
+/* the caller explicitly calls colamd_report or symamd_report. */
+#define PRINTF printf
+
+/* In C, matrices are 0-based and indices are reported as such in *_report */
+#define INDEX(i) (i)
+
+#endif /* MATLAB_MEX_FILE */
+
+/* ========================================================================== */
+/* === Prototypes of PRIVATE routines ======================================= */
+/* ========================================================================== */
+
+PRIVATE int init_rows_cols
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A [],
+    int p [],
+    int stats [COLAMD_STATS]
+) ;
+
+PRIVATE void init_scoring
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A [],
+    int head [],
+    double knobs [COLAMD_KNOBS],
+    int *p_n_row2,
+    int *p_n_col2,
+    int *p_max_deg
+) ;
+
+PRIVATE int find_ordering
+(
+    int n_row,
+    int n_col,
+    int Alen,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A [],
+    int head [],
+    int n_col2,
+    int max_deg,
+    int pfree
+) ;
+
+PRIVATE void order_children
+(
+    int n_col,
+    Colamd_Col Col [],
+    int p []
+) ;
+
+PRIVATE void detect_super_cols
+(
+
+#ifndef NDEBUG
+    int n_col,
+    Colamd_Row Row [],
+#endif /* NDEBUG */
+
+    Colamd_Col Col [],
+    int A [],
+    int head [],
+    int row_start,
+    int row_length
+) ;
+
+PRIVATE int garbage_collection
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A [],
+    int *pfree
+) ;
+
+PRIVATE int clear_mark
+(
+    int n_row,
+    Colamd_Row Row []
+) ;
+
+PRIVATE void print_report
+(
+    char *method,
+    int stats [COLAMD_STATS]
+) ;
+
+/* ========================================================================== */
+/* === Debugging prototypes and definitions ================================= */
+/* ========================================================================== */
+
+#ifndef NDEBUG
+
+/* colamd_debug is the *ONLY* global variable, and is only */
+/* present when debugging */
+
+PRIVATE int colamd_debug ;	/* debug print level */
+
+#define DEBUG0(params) { (void) PRINTF params ; }
+#define DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; }
+#define DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; }
+#define DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; }
+#define DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; }
+
+#ifdef MATLAB_MEX_FILE
+#define ASSERT(expression) (mxAssert ((expression), ""))
+#else
+#define ASSERT(expression) (assert (expression))
+#endif /* MATLAB_MEX_FILE */
+
+PRIVATE void colamd_get_debug	/* gets the debug print level from getenv */
+(
+    char *method
+) ;
+
+PRIVATE void debug_deg_lists
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int head [],
+    int min_score,
+    int should,
+    int max_deg
+) ;
+
+PRIVATE void debug_mark
+(
+    int n_row,
+    Colamd_Row Row [],
+    int tag_mark,
+    int max_mark
+) ;
+
+PRIVATE void debug_matrix
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A []
+) ;
+
+PRIVATE void debug_structures
+(
+    int n_row,
+    int n_col,
+    Colamd_Row Row [],
+    Colamd_Col Col [],
+    int A [],
+    int n_col2
+) ;
+
+#else /* NDEBUG */
+
+/* === No debugging ========================================================= */
+
+#define DEBUG0(params) ;
+#define DEBUG1(params) ;
+#define DEBUG2(params) ;
+#define DEBUG3(params) ;
+#define DEBUG4(params) ;
+
+#define ASSERT(expression) ((void) 0)
+
+#endif /* NDEBUG */
+
+/* ========================================================================== */
+
+
+
+/* ========================================================================== */
+/* === USER-CALLABLE ROUTINES: ============================================== */
+/* ========================================================================== */
+
+
+/* ========================================================================== */
+/* === colamd_recommended =================================================== */
+/* ========================================================================== */
+
+/*
+    The colamd_recommended routine returns the suggested size for Alen.  This
+    value has been determined to provide good balance between the number of
+    garbage collections and the memory requirements for colamd.  If any
+    argument is negative, a -1 is returned as an error condition.  This
+    function is also available as a macro defined in colamd.h, so that you
+    can use it for a statically-allocated array size.
+*/
+
+PUBLIC int colamd_recommended	/* returns recommended value of Alen. */
+(
+    /* === Parameters ======================================================= */
+
+    int nnz,			/* number of nonzeros in A */
+    int n_row,			/* number of rows in A */
+    int n_col			/* number of columns in A */
+)
+{
+    return (COLAMD_RECOMMENDED (nnz, n_row, n_col)) ; 
+}
+
+
+/* ========================================================================== */
+/* === colamd_set_defaults ================================================== */
+/* ========================================================================== */
+
+/*
+    The colamd_set_defaults routine sets the default values of the user-
+    controllable parameters for colamd:
+
+	knobs [0]	rows with knobs[0]*n_col entries or more are removed
+			prior to ordering in colamd.  Rows and columns with
+			knobs[0]*n_col entries or more are removed prior to
+			ordering in symamd and placed last in the output
+			ordering.
+
+	knobs [1]	columns with knobs[1]*n_row entries or more are removed
+			prior to ordering in colamd, and placed last in the
+			column permutation.  Symamd ignores this knob.
+
+	knobs [2..19]	unused, but future versions might use this
+*/
+
+PUBLIC void colamd_set_defaults
+(
+    /* === Parameters ======================================================= */
+
+    double knobs [COLAMD_KNOBS]		/* knob array */
+)
+{
+    /* === Local variables ================================================== */
+
+    int i ;
+
+    if (!knobs)
+    {
+	return ;			/* no knobs to initialize */
+    }
+    for (i = 0 ; i < COLAMD_KNOBS ; i++)
+    {
+	knobs [i] = 0 ;
+    }
+    knobs [COLAMD_DENSE_ROW] = 0.5 ;	/* ignore rows over 50% dense */
+    knobs [COLAMD_DENSE_COL] = 0.5 ;	/* ignore columns over 50% dense */
+}
+
+
+/* ========================================================================== */
+/* === symamd =============================================================== */
+/* ========================================================================== */
+
+PUBLIC int symamd			/* return TRUE if OK, FALSE otherwise */
+(
+    /* === Parameters ======================================================= */
+
+    int n,				/* number of rows and columns of A */
+    int A [],				/* row indices of A */
+    int p [],				/* column pointers of A */
+    int perm [],			/* output permutation, size n+1 */
+    double knobs [COLAMD_KNOBS],	/* parameters (uses defaults if NULL) */
+    int stats [COLAMD_STATS],		/* output statistics and error codes */
+    void * (*allocate) (size_t, size_t),
+    					/* pointer to calloc (ANSI C) or */
+					/* mxCalloc (for Matlab mexFunction) */
+    void (*release) (void *)
+    					/* pointer to free (ANSI C) or */
+    					/* mxFree (for Matlab mexFunction) */
+)
+{
+    /* === Local variables ================================================== */
+
+    int *count ;		/* length of each column of M, and col pointer*/
+    int *mark ;			/* mark array for finding duplicate entries */
+    int *M ;	      /* row indices of matrix M */
+    int Mlen ;			/* length of M */
+    int n_row ;			/* number of rows in M */
+    int nnz ;			/* number of entries in A */
+    int i ;			/* row index of A */
+    int j ;			/* column index of A */
+    int k ;			/* row index of M */ 
+    int mnz ;			/* number of nonzeros in M */
+    int pp ;			/* index into a column of A */
+    int last_row ;		/* last row seen in the current column */
+    int length ;		/* number of nonzeros in a column */
+
+    double cknobs [COLAMD_KNOBS] ;		/* knobs for colamd */
+    double default_knobs [COLAMD_KNOBS] ;	/* default knobs for colamd */
+    int cstats [COLAMD_STATS] ;			/* colamd stats */
+
+#ifndef NDEBUG
+    colamd_get_debug ("symamd") ;
+#endif /* NDEBUG */
+
+    /* === Check the input arguments ======================================== */
+
+    if (!stats)
+    {
+	DEBUG0 (("symamd: stats not present\n")) ;
+	return (FALSE) ;
+    }
+    for (i = 0 ; i < COLAMD_STATS ; i++)
+    {
+	stats [i] = 0 ;
+    }
+    stats [COLAMD_STATUS] = COLAMD_OK ;
+    stats [COLAMD_INFO1] = -1 ;
+    stats [COLAMD_INFO2] = -1 ;
+
+    if (!A)
+    {
+    	stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
+	DEBUG0 (("symamd: A not present\n")) ;
+	return (FALSE) ;
+    }
+
+    if (!p)		/* p is not present */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
+	DEBUG0 (("symamd: p not present\n")) ;
+    	return (FALSE) ;
+    }
+
+    if (n < 0)		/* n must be >= 0 */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
+	stats [COLAMD_INFO1] = n ;
+	DEBUG0 (("symamd: n negative %d\n", n)) ;
+    	return (FALSE) ;
+    }
+
+    nnz = p [n] ;
+    if (nnz < 0)	/* nnz must be >= 0 */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
+	stats [COLAMD_INFO1] = nnz ;
+	DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
+	return (FALSE) ;
+    }
+
+    if (p [0] != 0)
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
+	stats [COLAMD_INFO1] = p [0] ;
+	DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
+	return (FALSE) ;
+    }
+
+    /* === If no knobs, set default knobs =================================== */
+
+    if (!knobs)
+    {
+	colamd_set_defaults (default_knobs) ;
+	knobs = default_knobs ;
+    }
+
+    /* === Allocate count and mark ========================================== */
+
+    count = (int *) ((*allocate) (n+1, sizeof (int))) ;
+    if (!count)
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
+	DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
+	return (FALSE) ;
+    }
+
+    mark = (int *) ((*allocate) (n+1, sizeof (int))) ;
+    if (!mark)
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
+	(*release) ((void *) count) ;
+	DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
+	return (FALSE) ;
+    }
+
+    /* === Compute column counts of M, check if A is valid ================== */
+
+    stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
+
+    for (i = 0 ; i < n ; i++)
+    {
+    	mark [i] = -1 ;
+    }
+
+    for (j = 0 ; j < n ; j++)
+    {
+	last_row = -1 ;
+
+	length = p [j+1] - p [j] ;
+	if (length < 0)
+	{
+	    /* column pointers must be non-decreasing */
+	    stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
+	    stats [COLAMD_INFO1] = j ;
+	    stats [COLAMD_INFO2] = length ;
+	    (*release) ((void *) count) ;
+	    (*release) ((void *) mark) ;
+	    DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
+	    return (FALSE) ;
+	}
+
+	for (pp = p [j] ; pp < p [j+1] ; pp++)
+	{
+	    i = A [pp] ;
+	    if (i < 0 || i >= n)
+	    {
+		/* row index i, in column j, is out of bounds */
+		stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
+		stats [COLAMD_INFO1] = j ;
+		stats [COLAMD_INFO2] = i ;
+		stats [COLAMD_INFO3] = n ;
+		(*release) ((void *) count) ;
+		(*release) ((void *) mark) ;
+		DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
+		return (FALSE) ;
+	    }
+
+	    if (i <= last_row || mark [i] == j)
+	    {
+		/* row index is unsorted or repeated (or both), thus col */
+		/* is jumbled.  This is a notice, not an error condition. */
+		stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
+		stats [COLAMD_INFO1] = j ;
+		stats [COLAMD_INFO2] = i ;
+		(stats [COLAMD_INFO3]) ++ ;
+		DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
+	    }
+
+	    if (i > j && mark [i] != j)
+	    {
+		/* row k of M will contain column indices i and j */
+		count [i]++ ;
+		count [j]++ ;
+	    }
+
+	    /* mark the row as having been seen in this column */
+	    mark [i] = j ;
+
+	    last_row = i ;
+	}
+    }
+
+    if (stats [COLAMD_STATUS] == COLAMD_OK)
+    {
+	/* if there are no duplicate entries, then mark is no longer needed */
+	(*release) ((void *) mark) ;
+    }
+
+    /* === Compute column pointers of M ===================================== */
+
+    /* use output permutation, perm, for column pointers of M */
+    perm [0] = 0 ;
+    for (j = 1 ; j <= n ; j++)
+    {
+	perm [j] = perm [j-1] + count [j-1] ;
+    }
+    for (j = 0 ; j < n ; j++)
+    {
+	count [j] = perm [j] ;
+    }
+
+    /* === Construct M ====================================================== */
+
+    mnz = perm [n] ;
+    n_row = mnz / 2 ;
+    Mlen = colamd_recommended (mnz, n_row, n) ;
+    M = (int *) ((*allocate) (Mlen, sizeof (int))) ;
+    DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %d\n",
+    	n_row, n, mnz, Mlen)) ;
+
+    if (!M)
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
+	(*release) ((void *) count) ;
+	(*release) ((void *) mark) ;
+	DEBUG0 (("symamd: allocate M (size %d) failed\n", Mlen)) ;
+	return (FALSE) ;
+    }
+
+    k = 0 ;
+
+    if (stats [COLAMD_STATUS] == COLAMD_OK)
+    {
+	/* Matrix is OK */
+	for (j = 0 ; j < n ; j++)
+	{
+	    ASSERT (p [j+1] - p [j] >= 0) ;
+	    for (pp = p [j] ; pp < p [j+1] ; pp++)
+	    {
+		i = A [pp] ;
+		ASSERT (i >= 0 && i < n) ;
+		if (i > j)
+		{
+		    /* row k of M contains column indices i and j */
+		    M [count [i]++] = k ;
+		    M [count [j]++] = k ;
+		    k++ ;
+		}
+	    }
+	}
+    }
+    else
+    {
+	/* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
+	DEBUG0 (("symamd: Duplicates in A.\n")) ;
+	for (i = 0 ; i < n ; i++)
+	{
+	    mark [i] = -1 ;
+	}
+	for (j = 0 ; j < n ; j++)
+	{
+	    ASSERT (p [j+1] - p [j] >= 0) ;
+	    for (pp = p [j] ; pp < p [j+1] ; pp++)
+	    {
+		i = A [pp] ;
+		ASSERT (i >= 0 && i < n) ;
+		if (i > j && mark [i] != j)
+		{
+		    /* row k of M contains column indices i and j */
+		    M [count [i]++] = k ;
+		    M [count [j]++] = k ;
+		    k++ ;
+		    mark [i] = j ;
+		}
+	    }
+	}
+	(*release) ((void *) mark) ;
+    }
+
+    /* count and mark no longer needed */
+    (*release) ((void *) count) ;
+    ASSERT (k == n_row) ;
+
+    /* === Adjust the knobs for M =========================================== */
+
+    for (i = 0 ; i < COLAMD_KNOBS ; i++)
+    {
+	cknobs [i] = knobs [i] ;
+    }
+
+    /* there are no dense rows in M */
+    cknobs [COLAMD_DENSE_ROW] = 1.0 ;
+
+    if (n_row != 0 && n < n_row)
+    {
+	/* On input, the knob is a fraction of 1..n, the number of rows of A. */
+	/* Convert it to a fraction of 1..n_row, of the number of rows of M. */
+    	cknobs [COLAMD_DENSE_COL] = (knobs [COLAMD_DENSE_ROW] * n) / n_row ;
+    }
+    else
+    {
+	/* no dense columns in M */
+    	cknobs [COLAMD_DENSE_COL] = 1.0 ;
+    }
+
+    DEBUG0 (("symamd: dense col knob for M: %g\n", cknobs [COLAMD_DENSE_COL])) ;
+
+    /* === Order the columns of M =========================================== */
+
+    if (!colamd (n_row, n, Mlen, M, perm, cknobs, cstats))
+    {
+	/* This "cannot" happen, unless there is a bug in the code. */
+	stats [COLAMD_STATUS] = COLAMD_ERROR_internal_error ;
+	(*release) ((void *) M) ;
+	DEBUG0 (("symamd: internal error!\n")) ;
+	return (FALSE) ;
+    }
+
+    /* Note that the output permutation is now in perm */
+
+    /* === get the statistics for symamd from colamd ======================== */
+
+    /* note that a dense column in colamd means a dense row and col in symamd */
+    stats [COLAMD_DENSE_ROW]    = cstats [COLAMD_DENSE_COL] ;
+    stats [COLAMD_DENSE_COL]    = cstats [COLAMD_DENSE_COL] ;
+    stats [COLAMD_DEFRAG_COUNT] = cstats [COLAMD_DEFRAG_COUNT] ;
+
+    /* === Free M =========================================================== */
+
+    (*release) ((void *) M) ;
+    DEBUG0 (("symamd: done.\n")) ;
+    return (TRUE) ;
+
+}
+
+/* ========================================================================== */
+/* === colamd =============================================================== */
+/* ========================================================================== */
+
+/*
+    The colamd routine computes a column ordering Q of a sparse matrix
+    A such that the LU factorization P(AQ) = LU remains sparse, where P is
+    selected via partial pivoting.   The routine can also be viewed as
+    providing a permutation Q such that the Cholesky factorization
+    (AQ)'(AQ) = LL' remains sparse.
+*/
+
+PUBLIC int colamd		/* returns TRUE if successful, FALSE otherwise*/
+(
+    /* === Parameters ======================================================= */
+
+    int n_row,			/* number of rows in A */
+    int n_col,			/* number of columns in A */
+    int Alen,			/* length of A */
+    int A [],			/* row indices of A */
+    int p [],			/* pointers to columns in A */
+    double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
+    int stats [COLAMD_STATS]	/* output statistics and error codes */
+)
+{
+    /* === Local variables ================================================== */
+
+    int i ;			/* loop index */
+    int nnz ;			/* nonzeros in A */
+    int Row_size ;		/* size of Row [], in integers */
+    int Col_size ;		/* size of Col [], in integers */
+    int need ;			/* minimum required length of A */
+    Colamd_Row *Row ;		/* pointer into A of Row [0..n_row] array */
+    Colamd_Col *Col ;		/* pointer into A of Col [0..n_col] array */
+    int n_col2 ;		/* number of non-dense, non-empty columns */
+    int n_row2 ;		/* number of non-dense, non-empty rows */
+    int ngarbage ;		/* number of garbage collections performed */
+    int max_deg ;		/* maximum row degree */
+    double default_knobs [COLAMD_KNOBS] ;	/* default knobs array */
+
+#ifndef NDEBUG
+    colamd_get_debug ("colamd") ;
+#endif /* NDEBUG */
+
+    /* === Check the input arguments ======================================== */
+
+    if (!stats)
+    {
+	DEBUG0 (("colamd: stats not present\n")) ;
+	return (FALSE) ;
+    }
+    for (i = 0 ; i < COLAMD_STATS ; i++)
+    {
+	stats [i] = 0 ;
+    }
+    stats [COLAMD_STATUS] = COLAMD_OK ;
+    stats [COLAMD_INFO1] = -1 ;
+    stats [COLAMD_INFO2] = -1 ;
+
+    if (!A)		/* A is not present */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
+	DEBUG0 (("colamd: A not present\n")) ;
+	return (FALSE) ;
+    }
+
+    if (!p)		/* p is not present */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
+	DEBUG0 (("colamd: p not present\n")) ;
+    	return (FALSE) ;
+    }
+
+    if (n_row < 0)	/* n_row must be >= 0 */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
+	stats [COLAMD_INFO1] = n_row ;
+	DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
+    	return (FALSE) ;
+    }
+
+    if (n_col < 0)	/* n_col must be >= 0 */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
+	stats [COLAMD_INFO1] = n_col ;
+	DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
+    	return (FALSE) ;
+    }
+
+    nnz = p [n_col] ;
+    if (nnz < 0)	/* nnz must be >= 0 */
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
+	stats [COLAMD_INFO1] = nnz ;
+	DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
+	return (FALSE) ;
+    }
+
+    if (p [0] != 0)
+    {
+	stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero	;
+	stats [COLAMD_INFO1] = p [0] ;
+	DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
+	return (FALSE) ;
+    }
+
+    /* === If no knobs, set default knobs =================================== */
+
+    if (!knobs)
+    {
+	colamd_set_defaults (default_knobs) ;
+	knobs = default_knobs ;
+    }
+
+    /* === Allocate the Row and Col arrays from array A ===================== */
+
+    Col_size = COLAMD_C (n_col) ;
+    Row_size = COLAMD_R (n_row) ;
+    need = 2*nnz + n_col + Col_size + Row_size ;
+
+    if (need > Alen)
+    {
+	/* not enough space in array A to perform the ordering */
+	stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
+	stats [COLAMD_INFO1] = need ;
+	stats [COLAMD_INFO2] = Alen ;
+	DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
+	return (FALSE) ;
+    }
+
+    Alen -= Col_size + Row_size ;
+    Col = (Colamd_Col *) &A [Alen] ;
+    Row = (Colamd_Row *) &A [Alen + Col_size] ;
+
+    /* === Construct the row and column data structures ===================== */
+
+    if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
+    {
+	/* input matrix is invalid */
+	DEBUG0 (("colamd: Matrix invalid\n")) ;
+	return (FALSE) ;
+    }
+
+    /* === Initialize scores, kill dense rows/columns ======================= */
+
+    init_scoring (n_row, n_col, Row, Col, A, p, knobs,
+	&n_row2, &n_col2, &max_deg) ;
+
+    /* === Order the supercolumns =========================================== */
+
+    ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
+	n_col2, max_deg, 2*nnz) ;
+
+    /* === Order the non-principal columns ================================== */
+
+    order_children (n_col, Col, p) ;
+
+    /* === Return statistics in stats ======================================= */
+
+    stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
+    stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
+    stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
+    DEBUG0 (("colamd: done.\n")) ; 
+    return (TRUE) ;
+}
+
+
+/* ========================================================================== */
+/* === colamd_report ======================================================== */
+/* ========================================================================== */
+
+PUBLIC void colamd_report
+(
+    int stats [COLAMD_STATS]
+)
+{
+    print_report ("colamd", stats) ;
+}
+
+
+/* ========================================================================== */
+/* === symamd_report ======================================================== */
+/* ========================================================================== */
+
+PUBLIC void symamd_report
+(
+    int stats [COLAMD_STATS]
+)
+{
+    print_report ("symamd", stats) ;
+}
+
+
+
+/* ========================================================================== */
+/* === NON-USER-CALLABLE ROUTINES: ========================================== */
+/* ========================================================================== */
+
+/* There are no user-callable routines beyond this point in the file */
+
+
+/* ========================================================================== */
+/* === init_rows_cols ======================================================= */
+/* ========================================================================== */
+
+/*
+    Takes the column form of the matrix in A and creates the row form of the
+    matrix.  Also, row and column attributes are stored in the Col and Row
+    structs.  If the columns are un-sorted or contain duplicate row indices,
+    this routine will also sort and remove duplicate row indices from the
+    column form of the matrix.  Returns FALSE if the matrix is invalid,
+    TRUE otherwise.  Not user-callable.
+*/
+
+PRIVATE int init_rows_cols	/* returns TRUE if OK, or FALSE otherwise */
+(
+    /* === Parameters ======================================================= */
+
+    int n_row,			/* number of rows of A */
+    int n_col,			/* number of columns of A */
+    Colamd_Row Row [],		/* of size n_row+1 */
+    Colamd_Col Col [],		/* of size n_col+1 */
+    int A [],			/* row indices of A, of size Alen */
+    int p [],			/* pointers to columns in A, of size n_col+1 */
+    int stats [COLAMD_STATS]	/* colamd statistics */ 
+)
+{
+    /* === Local variables ================================================== */
+
+    int col ;			/* a column index */
+    int row ;			/* a row index */
+    int *cp ;			/* a column pointer */
+    int *cp_end ;		/* a pointer to the end of a column */
+    int *rp ;			/* a row pointer */
+    int *rp_end ;		/* a pointer to the end of a row */
+    int last_row ;		/* previous row */
+
+    /* === Initialize columns, and check column pointers ==================== */
+
+    for (col = 0 ; col < n_col ; col++)
+    {
+	Col [col].start = p [col] ;
+	Col [col].length = p [col+1] - p [col] ;
+
+	if (Col [col].length < 0)
+	{
+	    /* column pointers must be non-decreasing */
+	    stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
+	    stats [COLAMD_INFO1] = col ;
+	    stats [COLAMD_INFO2] = Col [col].length ;
+	    DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
+	    return (FALSE) ;
+	}
+
+	Col [col].shared1.thickness = 1 ;
+	Col [col].shared2.score = 0 ;
+	Col [col].shared3.prev = EMPTY ;
+	Col [col].shared4.degree_next = EMPTY ;
+    }
+
+    /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
+
+    /* === Scan columns, compute row degrees, and check row indices ========= */
+
+    stats [COLAMD_INFO3] = 0 ;	/* number of duplicate or unsorted row indices*/
+
+    for (row = 0 ; row < n_row ; row++)
+    {
+	Row [row].length = 0 ;
+	Row [row].shared2.mark = -1 ;
+    }
+
+    for (col = 0 ; col < n_col ; col++)
+    {
+	last_row = -1 ;
+
+	cp = &A [p [col]] ;
+	cp_end = &A [p [col+1]] ;
+
+	while (cp < cp_end)
+	{
+	    row = *cp++ ;
+
+	    /* make sure row indices within range */
+	    if (row < 0 || row >= n_row)
+	    {
+		stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
+		stats [COLAMD_INFO1] = col ;
+		stats [COLAMD_INFO2] = row ;
+		stats [COLAMD_INFO3] = n_row ;
+		DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
+		return (FALSE) ;
+	    }
+
+	    if (row <= last_row || Row [row].shared2.mark == col)
+	    {
+		/* row index are unsorted or repeated (or both), thus col */
+		/* is jumbled.  This is a notice, not an error condition. */
+		stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
+		stats [COLAMD_INFO1] = col ;
+		stats [COLAMD_INFO2] = row ;
+		(stats [COLAMD_INFO3]) ++ ;
+		DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
+	    }
+
+	    if (Row [row].shared2.mark != col)
+	    {
+		Row [row].length++ ;
+	    }
+	    else
+	    {
+		/* this is a repeated entry in the column, */
+		/* it will be removed */
+		Col [col].length-- ;
+	    }
+
+	    /* mark the row as having been seen in this column */
+	    Row [row].shared2.mark = col ;
+
+	    last_row = row ;
+	}
+    }
+
+    /* === Compute row pointers ============================================= */
+
+    /* row form of the matrix starts directly after the column */
+    /* form of matrix in A */
+    Row [0].start = p [n_col] ;
+    Row [0].shared1.p = Row [0].start ;
+    Row [0].shared2.mark = -1 ;
+    for (row = 1 ; row < n_row ; row++)
+    {
+	Row [row].start = Row [row-1].start + Row [row-1].length ;
+	Row [row].shared1.p = Row [row].start ;
+	Row [row].shared2.mark = -1 ;
+    }
+
+    /* === Create row form ================================================== */
+
+    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
+    {
+	/* if cols jumbled, watch for repeated row indices */
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    cp = &A [p [col]] ;
+	    cp_end = &A [p [col+1]] ;
+	    while (cp < cp_end)
+	    {
+		row = *cp++ ;
+		if (Row [row].shared2.mark != col)
+		{
+		    A [(Row [row].shared1.p)++] = col ;
+		    Row [row].shared2.mark = col ;
+		}
+	    }
+	}
+    }
+    else
+    {
+	/* if cols not jumbled, we don't need the mark (this is faster) */
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    cp = &A [p [col]] ;
+	    cp_end = &A [p [col+1]] ;
+	    while (cp < cp_end)
+	    {
+		A [(Row [*cp++].shared1.p)++] = col ;
+	    }
+	}
+    }
+
+    /* === Clear the row marks and set row degrees ========================== */
+
+    for (row = 0 ; row < n_row ; row++)
+    {
+	Row [row].shared2.mark = 0 ;
+	Row [row].shared1.degree = Row [row].length ;
+    }
+
+    /* === See if we need to re-create columns ============================== */
+
+    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
+    {
+    	DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
+
+#ifndef NDEBUG
+	/* make sure column lengths are correct */
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    p [col] = Col [col].length ;
+	}
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    rp = &A [Row [row].start] ;
+	    rp_end = rp + Row [row].length ;
+	    while (rp < rp_end)
+	    {
+		p [*rp++]-- ;
+	    }
+	}
+	for (col = 0 ; col < n_col ; col++)
+	{
+	    ASSERT (p [col] == 0) ;
+	}
+	/* now p is all zero (different than when debugging is turned off) */
+#endif /* NDEBUG */
+
+	/* === Compute col pointers ========================================= */
+
+	/* col form of the matrix starts at A [0]. */
+	/* Note, we may have a gap between the col form and the row */
+	/* form if there were duplicate entries, if so, it will be */
+	/* removed upon the first garbage collection */
+	Col [0].start = 0 ;
+	p [0] = Col [0].start ;
+	for (col = 1 ; col < n_col ; col++)
+	{
+	    /* note that the lengths here are for pruned columns, i.e. */
+	    /* no duplicate row indices will exist for these columns */
+	    Col [col].start = Col [col-1].start + Col [col-1].length ;
+	    p [col] = Col [col].start ;
+	}
+
+	/* === Re-create col form =========================================== */
+
+	for (row = 0 ; row < n_row ; row++)
+	{
+	    rp = &A [Row [row].start] ;
+	    rp_end = rp + Row [row].length ;
+	    while (rp < rp_end)
+	    {
+		A [(p [*rp++])++] = row ;
+	    }
+	}
+    }
+
+    /* === Done.  Matrix is not (or no longer) jumbled ====================== */
+
+    return (TRUE) ;
+}
+
+
+/* ========================================================================== */
+/* === init_scoring ========================================================= */
+/* ========================================================================== */
+
+/*
+    Kills dense or empty columns and rows, calculates an initial score for
+    each column, and places all columns in the degree lists.  Not user-callable.
+*/
+
+PRIVATE void init_scoring
+(
+    /* === Parameters ======================================================= */
+
+    int n_row,			/* number of rows of A */
+    int n_col,			/* number of columns of A */
+    Colamd_Row Row [],		/* of size n_row+1 */
+    Colamd_Col Col [],		/* of size n_col+1 */
+    int A [],			/* column form and row form of A */
+    int head [],		/* of size n_col+1 */
+    double knobs [COLAMD_KNOBS],/* parameters */
+    int *p_n_row2,		/* number of non-dense, non-empty rows */
+    int *p_n_col2,		/* number of non-dense, non-empty columns */
+    int *p_max_deg		/* maximum row degree */
+)
+{
+    /* === Local variables ================================================== */
+
+    int c ;         /* a column index */
+    int r, row ;		/* a row index */
+    int *cp ;		    /* a column pointer */
+    int deg ;	      /* degree of a row or column */
+    int *cp_end ;	  /* a pointer to the end of a column */
+    int *new_cp ;	  /* new column pointer */
+    int col_length ;		/* length of pruned column */
+    int score ;	    /* current column score */
+    int n_col2 ;    /* number of non-dense, non-empty columns */
+    int n_row2 ;    /* number of non-dense, non-empty rows */
+    int dense_row_count ;	/* remove rows with more entries than this */
+    int dense_col_count ;	/* remove cols with more entries than this */
+    int min_score ; /* smallest column score */
+    int max_deg ;	  /* maximum row degree */
+    int next_col ;  /* Used to add to degree list.*/
+
+#ifndef NDEBUG
+    int debug_count ;		/* debug only. */
+#endif /* NDEBUG */
+
+    /* === Extract knobs ==================================================== */
+
+    dense_row_count = (int) MAX (0, MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ;
+    dense_col_count = (int) MAX (0, MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ;
+    DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
+    max_deg = 0 ;
+    n_col2 = n_col ;
+    n_row2 = n_row ;
+
+    /* === Kill empty columns =============================================== */
+
+    /* Put the empty columns at the end in their natural order, so that LU */
+    /* factorization can proceed as far as possible. */
+    for (c = n_col-1 ; c >= 0 ; c--)
+    {
+	deg = Col [c].length ;
+	if (deg == 0)
+	{
+	    /* this is a empty column, kill and order it last */
+	    Col [c].shared2.order = --n_col2 ;
+	    KILL_PRINCIPAL_COL (c) ;
+	}
+    }
+    DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
+
+    /* === Kill dense columns =============================================== */
+
+    /* Put the dense columns at the end, in their natural order */
+    for (c = n_col-1 ; c >= 0 ; c--)
+    {
+	/* skip any dead columns */
+	if (COL_IS_DEAD (c))
+	{
+	    continue ;
+	}
+	deg = Col [c].length ;
+	if (deg > dense_col_count)
+	{
+	    /* this is a dense column, kill and order it last */
+	    Col [c].shared2.order = --n_col2 ;
+	    /* decrement the row degrees */
+	    cp = &A [Col [c].start] ;
+	    cp_end = cp + Col [c].length ;
+	    while (cp < cp_end)
+	    {
+		Row [*cp++].shared1.degree-- ;
+	    }
+	    KILL_PRINCIPAL_COL (c) ;
+	}
+    }
+    DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
+
+    /* === Kill dense and empty rows ======================================== */
+
+    for (r = 0 ; r < n_row ; r++)
+    {
+	deg = Row [r].shared1.degree ;
+	ASSERT (deg >= 0 && deg <= n_col) ;
+	if (deg > dense_row_count || deg == 0)
+	{
+	    /* kill a dense or empty row */
+	    KILL_ROW (r) ;
+	    --n_row2 ;
+	}
+	else
+	{
+	    /* keep track of max degree of remaining rows */
+	    max_deg = MAX (max_deg, deg) ;
+	}
+    }
+    DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
+
+    /* === Compute initial column scores ==================================== */
+
+    /* At this point the row degrees are accurate.  They reflect the number */
+    /* of "live" (non-dense) columns in each row.  No empty rows exist. */
+    /* Some "live" columns may contain only dead rows, however.  These are */
+    /* pruned in the code below. */
+
+    /* now find the initial matlab score for each column */
+    for (c = n_col-1 ; c >= 0 ; c--)
+    {
+	/* skip dead column */
+	if (COL_IS_DEAD (c))
+	{
+	    continue ;
+	}
+	score = 0 ;
+	cp = &A [Col [c].start] ;
+	new_cp = cp ;
+	cp_end = cp + Col [c].length ;
+	while (cp < cp_end)
+	{
+	    /* get a row */
+	    row = *cp++ ;
+	    /* skip if dead */
+	    if (ROW_IS_DEAD (row))
+	    {
+		continue ;
+	    }
+	    /* compact the column */
+	    *new_cp++ = row ;
+	    /* add row's external degree */
+	    score += Row [row].shared1.degree - 1 ;
+	    /* guard against integer overflow */
+	    score = MIN (score, n_col) ;
+	}
+	/* determine pruned column length */
+	col_length = (int) (new_cp - &A [Col [c].start]) ;
+	if (col_length == 0)
+	{
+	    /* a newly-made null column (all rows in this col are "dense" */
+	    /* and have already been killed) */
+	    DEBUG2 (("Newly null killed: %d\n", c)) ;
+	    Col [c].shared2.order = --n_col2 ;
+	    KILL_PRINCIPAL_COL (c) ;
+	}
+	else
+	{
+	    /* set column length and set score */
+	    ASSERT (score >= 0) ;
+	    ASSERT (score <= n_col) ;
+	    Col [c].length = col_length ;
+	    Col [c].shared2.score = score ;
+	}
+    }
+    DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
+    	n_col-n_col2)) ;
+
+    /* At this point, all empty rows and columns are dead.  All live columns */
+    /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
+    /* yet).  Rows may contain dead columns, but all live rows contain at */
+    /* least one live column. */
+
+#ifndef NDEBUG
+    debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
+#endif /* NDEBUG */
+
+    /* === Initialize degree lists ========================================== */
+
+#ifndef NDEBUG
+    debug_count = 0 ;
+#endif /* NDEBUG */
+
+    /* clear the hash buckets */
+    for (c = 0 ; c <= n_col ; c++)
+    {
+	head [c] = EMPTY ;
+    }
+    min_score = n_col ;
+    /* place in reverse order, so low column indices are at the front */
+    /* of the lists.  This is to encourage natural tie-breaking */
+    for (c = n_col-1 ; c >= 0 ; c--)
+    {
+	/* only add principal columns to degree lists */
+	if (COL_IS_ALIVE (c))
+	{
+	    DEBUG4 (("place %d score %d minscore %d ncol %d\n",
+		c, Col [c].shared2.score, min_score, n_col)) ;
+
+	    /* === Add columns score to DList =============================== */
+
+	    score = Col [c].shared2.score ;
+
+	    ASSERT (min_score >= 0) ;
+	    ASSERT (min_score <= n_col) ;
+	    ASSERT (score >= 0) ;
+	    ASSERT (score <= n_col) ;
+	    ASSERT (head [score] >= EMPTY) ;
+
+	    /* now add this column to dList at proper score location */
+	    next_col = head [score] ;
+	    Col [c].shared3.prev = EMPTY ;
+	    Col [c].shared4.degree_next = next_col ;
+
+	    /* if there already was a column with the same score, set its */
+	    /* previous pointer to this new column */
+	    if (next_col != EMPTY)
+	    {
+		Col [next_col].shared3.prev = c ;
+	    }
+	    head [score] = c ;
+
+	    /* see if this score is less than current min */
+	    min_score = MIN (min_score, score) ;
+
+#ifndef NDEBUG
+	    debug_count++ ;
+#endif /* NDEBUG */
+
+	}
+    }
+
+#ifndef NDEBUG
+    DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
+	debug_count, n_col, n_col-debug_count)) ;
+    ASSERT (debug_count == n_col2) ;
+    debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
+#endif /* NDEBUG */
+
+    /* === Return number of remaining columns, and max row degree =========== */
+
+    *p_n_col2 = n_col2 ;
+    *p_n_row2 = n_row2 ;
+    *p_max_deg = max_deg ;
+}
+
+
+/* ========================================================================== */
+/* === find_ordering ======================================================== */
+/* ========================================================================== */
+
+/*
+    Order the principal columns of the supercolumn form of the matrix
+    (no supercolumns on input).  Uses a minimum approximate column minimum
+    degree ordering method.  Not user-callable.
+*/
+
+PRIVATE int find_ordering	/* return the number of garbage collections */
+(
+    /* === Parameters ======================================================= */
+
+    int n_row,			/* number of rows of A */
+    int n_col,			/* number of columns of A */
+    int Alen,			/* size of A, 2*nnz + n_col or larger */
+    Colamd_Row Row [],		/* of size n_row+1 */
+    Colamd_Col Col [],		/* of size n_col+1 */
+    int A [],			/* column form and row form of A */
+    int head [],		/* of size n_col+1 */
+    int n_col2,			/* Remaining columns to order */
+    int max_deg,		/* Maximum row degree */
+    int pfree			/* index of first free slot (2*nnz on entry) */
+)
+{
+    /* === Local variables ================================================== */
+
+    int k ;			/* current pivot ordering step */
+    int pivot_col ;		/* current pivot column */
+    int *cp ;			/* a column pointer */
+    int *rp ;			/* a row pointer */
+    int pivot_row ;		/* current pivot row */
+    int *new_cp ;		/* modified column pointer */
+    int *new_rp ;		/* modified row pointer */
+    int pivot_row_start ;	/* pointer to start of pivot row */
+    int pivot_row_degree ;	/* number of columns in pivot row */
+    int pivot_row_length ;	/* number of supercolumns in pivot row */
+    int pivot_col_score ;	/* score of pivot column */
+    int needed_memory ;		/* free space needed for pivot row */
+    int *cp_end ;		/* pointer to the end of a column */
+    int *rp_end ;		/* pointer to the end of a row */
+    int row ;			/* a row index */
+    int col ;			/* a column index */
+    int max_score ;		/* maximum possible score */
+    int cur_score ;		/* score of current column */
+    unsigned int hash ;		/* hash value for supernode detection */
+    int head_column ;		/* head of hash bucket */
+    int first_col ;		/* first column in hash bucket */
+    int tag_mark ;		/* marker value for mark array */
+    int row_mark ;		/* Row [row].shared2.mark */
+    int set_difference ;	/* set difference size of row with pivot row */
+    int min_score ;		/* smallest column score */
+    int col_thickness ;		/* "thickness" (no. of columns in a supercol) */
+    int max_mark ;		/* maximum value of tag_mark */
+    int pivot_col_thickness ;	/* number of columns represented by pivot col */
+    int prev_col ;		/* Used by Dlist operations. */
+    int next_col ;		/* Used by Dlist operations. */
+    int ngarbage ;		/* number of garbage collections performed */
+
+#ifndef NDEBUG
+    int debug_d ;		/* debug loop counter */
+    int debug_step = 0 ;	/* debug loop counter */
+#endif /* NDEBUG */
+
+    /* === Initialization and clear mark ==================================== */
+
+    max_mark = INT_MAX - n_col ;	/* INT_MAX defined in <limits.h> */
+    tag_mark = clear_mark (n_row, Row) ;
+    min_score = 0 ;
+    ngarbage = 0 ;
+    DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
+
+    /* === Order the columns ================================================ */
+
+    for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
+    {
+
+#ifndef NDEBUG
+	if (debug_step % 100 == 0)
+	{
+	    DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
+	}
+	else
+	{
+	    DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
+	}
+	debug_step++ ;
+	debug_deg_lists (n_row, n_col, Row, Col, head,
+		min_score, n_col2-k, max_deg) ;
+	debug_matrix (n_row, n_col, Row, Col, A) ;
+#endif /* NDEBUG */
+
+	/* === Select pivot column, and order it ============================ */
+
+	/* make sure degree list isn't empty */
+	ASSERT (min_score >= 0) ;
+	ASSERT (min_score <= n_col) ;
+	ASSERT (head [min_score] >= EMPTY) ;
+
+#ifndef NDEBUG
+	for (debug_d = 0 ; debug_d < min_score ; debug_d++)
+	{
+	    ASSERT (head [debug_d] == EMPTY) ;
+	}
+#endif /* NDEBUG */
+
+	/* get pivot column from head of minimum degree list */
+	while (head [min_score] == EMPTY && min_score < n_col)
+	{
+	    min_score++ ;
+	}
+	pivot_col = head [min_score] ;
+	ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
+	next_col = Col [pivot_col].shared4.degree_next ;
+	head [min_score] = next_col ;
+	if (next_col != EMPTY)
+	{
+	    Col [next_col].shared3.prev = EMPTY ;
+	}
+
+	ASSERT (COL_IS_ALIVE (pivot_col)) ;
+	DEBUG3 (("Pivot col: %d\n", pivot_col)) ;
+
+	/* remember score for defrag check */
+	pivot_col_score = Col [pivot_col].shared2.score ;
+
+	/* the pivot column is the kth column in the pivot order */
+	Col [pivot_col].shared2.order = k ;
+
+	/* increment order count by column thickness */
+	pivot_col_thickness = Col [pivot_col].shared1.thickness ;
+	k += pivot_col_thickness ;
+	ASSERT (pivot_col_thickness > 0) ;
+
+	/* === Garbage_collection, if necessary ============================= */
+
+	needed_memory = MIN (pivot_col_score, n_col - k) ;
+	if (pfree + needed_memory >= Alen)
+	{
+	    pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
+	    ngarbage++ ;
+	    /* after garbage collection we will have enough */
+	    ASSERT (pfree + needed_memory < Alen) ;
+	    /* garbage collection has wiped out the Row[].shared2.mark array */
+	    tag_mark = clear_mark (n_row, Row) ;
+
+#ifndef NDEBUG
+	    debug_matrix (n_row, n_col, Row, Col, A) ;
+#endif /* NDEBUG */
+	}
+
+	/* === Compute pivot row pattern ==================================== */
+
+	/* get starting location for this new merged row */
+	pivot_row_start = pfree ;
+
+	/* initialize new row counts to zero */
+	pivot_row_degree = 0 ;
+
+	/* tag pivot column as having been visited so it isn't included */
+	/* in merged pivot row */
+	Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
+
+	/* pivot row is the union of all rows in the pivot column pattern */
+	cp = &A [Col [pivot_col].start] ;
+	cp_end = cp + Col [pivot_col].length ;
+	while (cp < cp_end)
+	{
+	    /* get a row */
+	    row = *cp++ ;
+	    DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
+	    /* skip if row is dead */
+	    if (ROW_IS_DEAD (row))
+	    {
+		continue ;
+	    }
+	    rp = &A [Row [row].start] ;
+	    rp_end = rp + Row [row].length ;
+	    while (rp < rp_end)
+	    {
+		/* get a column */
+		col = *rp++ ;
+		/* add the column, if alive and untagged */
+		col_thickness = Col [col].shared1.thickness ;
+		if (col_thickness > 0 && COL_IS_ALIVE (col))
+		{
+		    /* tag column in pivot row */
+		    Col [col].shared1.thickness = -col_thickness ;
+		    ASSERT (pfree < Alen) ;
+		    /* place column in pivot row */
+		    A [pfree++] = col ;
+		    pivot_row_degree += col_thickness ;
+		}