/* ========================================================================== */

 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */

 /* ========================================================================== */

 /* COLAMD / SYMAMD

     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.

         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 at cise.ufl.edu), University of Florida.  The algorithm was

         developed in collaboration with John Gilbert, Xerox PARC, and Esmond

         Ng, Oak Ridge National Laboratory.

     Acknowledgements:

         This work was supported by the National Science Foundation, under

         grants DMS-9504974 and DMS-9803599.

     Copyright and License:

         Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.

         COLAMD is also available under alternate licenses, contact T. Davis

         for details.

         This library is free software; you can redistribute it and/or

         modify it under the terms of the GNU Lesser General Public

         License as published by the Free Software Foundation; either

         version 2.1 of the License, or (at your option) any later version.

         This library is distributed in the hope that it will be useful,

         but WITHOUT ANY WARRANTY; without even the implied warranty of

         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

         Lesser General Public License for more details.

         You should have received a copy of the GNU Lesser General Public

         License along with this library; if not, write to the Free Software

         Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301

         USA

         Permission is hereby granted to use or copy this program under the

         terms of the GNU LGPL, provided that the Copyright, this License,

         and the Availability of the original version is retained on all copies.

         User documentation of any code that uses this code or any modified

         version of this code must cite the Copyright, this License, the

         Availability note, and "Used by permission." Permission to modify

         the code and to distribute modified code is granted, provided the

         Copyright, this License, and the Availability note are retained,

         and a notice that the code was modified is included.

     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.

         Appears as ACM Algorithm 836.

     See the ChangeLog file for changes since Version 1.0.

     References:

         T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column

         minimum degree ordering algorithm, ACM Transactions on Mathematical

         Software, vol. 30, no. 3., pp. 353-376, 2004.

         T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,

         an approximate column minimum degree ordering algorithm, ACM

         Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,

         2004.

 */

 /* ========================================================================== */

 /* === Description of user-callable routines ================================ */

 /* ========================================================================== */

 /* COLAMD includes both int and UF_long versions of all its routines.  The

  * description below is for the int version.  For UF_long, all int arguments

  * become UF_long.  UF_long is normally defined as long, except for WIN64.

     ----------------------------------------------------------------------------

     colamd_recommended:

     ----------------------------------------------------------------------------

         C syntax:

             #include "colamd.h"

             size_t colamd_recommended (int nnz, int n_row, int n_col) ;

             size_t colamd_l_recommended (UF_long nnz, UF_long n_row,

                 UF_long n_col) ;

         Purpose:

             Returns recommended value of Alen for use by colamd.  Returns 0

             if any input argument is negative.  The use of this routine

             is optional.  Not needed for symamd, which dynamically allocates

             its own memory.

             Note that in v2.4 and earlier, these routines returned int or long.

             They now return a value of type size_t.

         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 (double knobs [COLAMD_KNOBS]) ;

             colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;

         Purpose:

             Sets the default parameters.  The use of this routine is optional.

         Arguments:

             double knobs [COLAMD_KNOBS] ;       Output only.

                 NOTE: the meaning of the dense row/col knobs has changed in v2.4

                 knobs [0] and knobs [1] control dense row and col detection:

                 Colamd: rows with more than

                 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))

                 entries are removed prior to ordering.  Columns with more than

                 max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))

                 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

                 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (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 10.  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.

             knobs [2]: aggressive absorption

                 knobs [COLAMD_AGGRESSIVE] controls whether or not to do

                 aggressive absorption during the ordering.  Default is TRUE.

     ----------------------------------------------------------------------------

     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]) ;

             UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,

                 UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],

                 UF_long 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)

                 will be sufficient.  Note: the macro version does not check

                 for integer overflow, and thus is not recommended.  Use

                 the colamd_recommended routine instead.

             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 100

                 int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;

                 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,

                 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],

                 void (*allocate) (size_t, size_t), void (*release) (void *)) ;

             UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,

                 double knobs [COLAMD_KNOBS], UF_long 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).

                 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]) ;

             colamd_l_report (UF_long 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]) ;

             symamd_l_report (UF_long 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

 /* turn on debugging by uncommenting the following line

  #undef 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...

     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"

 #if 0 /* by mao */

 #include <limits.h>

 #include <math.h>

 #ifdef MATLAB_MEX_FILE

 #include "mex.h"

 #include "matrix.h"

 #endif /* MATLAB_MEX_FILE */

 #if !defined (NPRINT) || !defined (NDEBUG)

 #include <stdio.h>

 #endif

 #ifndef NULL

 #define NULL ((void *) 0)

 #endif

 #endif

 /* ========================================================================== */

 /* === int or UF_long ======================================================= */

 /* ========================================================================== */

 #if 0 /* by mao */

 /* define UF_long */

 #include "UFconfig.h"

 #endif

 #ifdef DLONG

 #define Int UF_long

 #define ID  UF_long_id

 #define Int_MAX UF_long_max

 #define COLAMD_recommended colamd_l_recommended

 #define COLAMD_set_defaults colamd_l_set_defaults

 #define COLAMD_MAIN colamd_l

 #define SYMAMD_MAIN symamd_l

 #define COLAMD_report colamd_l_report

 #define SYMAMD_report symamd_l_report

 #else

 #define Int int

 #define ID "%d"

 #define Int_MAX INT_MAX

 #define COLAMD_recommended colamd_recommended

 #define COLAMD_set_defaults colamd_set_defaults

 #define COLAMD_MAIN colamd

 #define SYMAMD_MAIN symamd

 #define COLAMD_report colamd_report

 #define SYMAMD_report symamd_report

 #endif

 /* ========================================================================== */

 /* === Row and Column structures ============================================ */

 /* ========================================================================== */

 /* User code that makes use of the colamd/symamd routines need not directly */

 /* reference these structures.  They are used only for colamd_recommended. */

 typedef struct Colamd_Col_struct

 {

     Int start ;         /* index for A of first row in this column, or DEAD */

                         /* if column is dead */

     Int length ;        /* number of rows in this column */

     union

     {

         Int thickness ; /* number of original columns represented by this */

                         /* col, if the column is alive */

         Int parent ;    /* parent in parent tree super-column structure, if */

                         /* the column is dead */

     } shared1 ;

     union

     {

         Int score ;     /* the score used to maintain heap, if col is alive */

         Int order ;     /* pivot ordering of this column, if col is dead */

     } shared2 ;

     union

     {

         Int headhash ;  /* head of a hash bucket, if col is at the head of */

                         /* a degree list */

         Int hash ;      /* hash value, if col is not in a degree list */

         Int prev ;      /* previous column in degree list, if col is in a */

                         /* degree list (but not at the head of a degree list) */

     } shared3 ;

     union

     {

         Int degree_next ;       /* next column, if col is in a degree list */

         Int hash_next ;         /* next column, if col is in a hash list */

     } shared4 ;

 } Colamd_Col ;

 typedef struct Colamd_Row_struct

 {

     Int start ;         /* index for A of first col in this row */

     Int length ;        /* number of principal columns in this row */

     union

     {

         Int degree ;    /* number of principal & non-principal columns in row */

         Int p ;         /* used as a row pointer in init_rows_cols () */

     } shared1 ;

     union

     {

         Int mark ;      /* for computing set differences and marking dead rows*/

         Int first_column ;/* first column in row (used in garbage collection) */

     } shared2 ;

 } Colamd_Row ;

 /* ========================================================================== */

 /* === Definitions ========================================================== */

 /* ========================================================================== */

 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */

 #define PUBLIC

 #define PRIVATE static

 #define DENSE_DEGREE(alpha,n) \

     ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))

 #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 =========================================== */

 /* ========================================================================== */

 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)

 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */

 #define INDEX(i) ((i)+1)

 #else

 /* In C, matrices are 0-based and indices are reported as such in *_report */

 #define INDEX(i) (i)

 #endif

 /* All output goes through the PRINTF macro.  */

 #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }

 /* ========================================================================== */

 /* === 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,

     Int aggressive

 ) ;

 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 tag_mark,

     Int max_mark,

     Int n_row,

     Colamd_Row Row []

 ) ;

 PRIVATE void print_report

 (

     char *method,

     Int stats [COLAMD_STATS]

 ) ;

 /* ========================================================================== */

 /* === Debugging prototypes and definitions ================================= */

 /* ========================================================================== */

 #ifndef NDEBUG

 #if 0 /* by mao */

 #include <assert.h>

 #endif

 /* colamd_debug is the *ONLY* global variable, and is only */

 /* present when debugging */

 PRIVATE Int colamd_debug = 0 ;  /* debug print level */

 #define DEBUG0(params) { PRINTF (params) ; }

 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }

 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }

 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }

 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }

 #if 0 /* by mao */

 #ifdef MATLAB_MEX_FILE

 #define ASSERT(expression) (mxAssert ((expression), ""))

 #else

 #define ASSERT(expression) (assert (expression))

 #endif /* MATLAB_MEX_FILE */

 #else

 #define ASSERT xassert

 #endif

 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)

 #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, or if integer overflow occurs, a 0 is returned as an

     error condition.  2*nnz space is required for the row and column

     indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is

     required for the Col and Row arrays, respectively, which are internal to

     colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the

     minimal amount of "elbow room", and nnz/5 more space is recommended for

     run time efficiency.

     Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.

     This function is not needed when using symamd.

 */

 /* add two values of type size_t, and check for integer overflow */

 static size_t t_add (size_t a, size_t b, int *ok)

 {

     (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;

     return ((*ok) ? (a + b) : 0) ;

 }

 /* compute a*k where k is a small integer, and check for integer overflow */

 static size_t t_mult (size_t a, size_t k, int *ok)

 {

     size_t i, s = 0 ;

     for (i = 0 ; i < k ; i++)

     {

         s = t_add (s, a, ok) ;

     }

     return (s) ;

 }

 /* size of the Col and Row structures */

 #define COLAMD_C(n_col,ok) \

     ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))

 #define COLAMD_R(n_row,ok) \

     ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))

 PUBLIC size_t 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 */

 )

 {

     size_t s, c, r ;

     int ok = TRUE ;

     if (nnz < 0 || n_row < 0 || n_col < 0)

     {

         return (0) ;

     }

     s = t_mult (nnz, 2, &ok) ;      /* 2*nnz */

     c = COLAMD_C (n_col, &ok) ;     /* size of column structures */

     r = COLAMD_R (n_row, &ok) ;     /* size of row structures */

     s = t_add (s, c, &ok) ;

     s = t_add (s, r, &ok) ;

     s = t_add (s, n_col, &ok) ;     /* elbow room */

     s = t_add (s, nnz/5, &ok) ;     /* elbow room */

     ok = ok && (s < Int_MAX) ;

     return (ok ? s : 0) ;

 }

 /* ========================================================================== */

 /* === colamd_set_defaults ================================================== */

 /* ========================================================================== */

 /*

     The colamd_set_defaults routine sets the default values of the user-

     controllable parameters for colamd and symamd:

         Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))

         entries are removed prior to ordering.  Columns with more than

         max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed

         prior to ordering, and placed last in the output column ordering. 

         Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))

         entries are removed prior to ordering, and placed last in the

         output ordering.

         knobs [0]       dense row control

         knobs [1]       dense column control

         knobs [2]       if nonzero, do aggresive absorption

         knobs [3..19]   unused, but future versions might use this

 */

 PUBLIC void COLAMD_set_defaults

 (

     /* === Parameters ======================================================= */

     double knobs [COLAMD_KNOBS]         /* knob array */

 )

 {