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

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