1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/colamd/colamd.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,3622 @@
     1.4 +/* ========================================================================== */
     1.5 +/* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
     1.6 +/* ========================================================================== */
     1.7 +
     1.8 +/* COLAMD / SYMAMD
     1.9 +
    1.10 +    colamd:  an approximate minimum degree column ordering algorithm,
    1.11 +        for LU factorization of symmetric or unsymmetric matrices,
    1.12 +        QR factorization, least squares, interior point methods for
    1.13 +        linear programming problems, and other related problems.
    1.14 +
    1.15 +    symamd:  an approximate minimum degree ordering algorithm for Cholesky
    1.16 +        factorization of symmetric matrices.
    1.17 +
    1.18 +    Purpose:
    1.19 +
    1.20 +        Colamd computes a permutation Q such that the Cholesky factorization of
    1.21 +        (AQ)'(AQ) has less fill-in and requires fewer floating point operations
    1.22 +        than A'A.  This also provides a good ordering for sparse partial
    1.23 +        pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
    1.24 +        factorization, and P is computed during numerical factorization via
    1.25 +        conventional partial pivoting with row interchanges.  Colamd is the
    1.26 +        column ordering method used in SuperLU, part of the ScaLAPACK library.
    1.27 +        It is also available as built-in function in MATLAB Version 6,
    1.28 +        available from MathWorks, Inc. (http://www.mathworks.com).  This
    1.29 +        routine can be used in place of colmmd in MATLAB.
    1.30 +
    1.31 +        Symamd computes a permutation P of a symmetric matrix A such that the
    1.32 +        Cholesky factorization of PAP' has less fill-in and requires fewer
    1.33 +        floating point operations than A.  Symamd constructs a matrix M such
    1.34 +        that M'M has the same nonzero pattern of A, and then orders the columns
    1.35 +        of M using colmmd.  The column ordering of M is then returned as the
    1.36 +        row and column ordering P of A. 
    1.37 +
    1.38 +    Authors:
    1.39 +
    1.40 +        The authors of the code itself are Stefan I. Larimore and Timothy A.
    1.41 +        Davis (davis at cise.ufl.edu), University of Florida.  The algorithm was
    1.42 +        developed in collaboration with John Gilbert, Xerox PARC, and Esmond
    1.43 +        Ng, Oak Ridge National Laboratory.
    1.44 +
    1.45 +    Acknowledgements:
    1.46 +
    1.47 +        This work was supported by the National Science Foundation, under
    1.48 +        grants DMS-9504974 and DMS-9803599.
    1.49 +
    1.50 +    Copyright and License:
    1.51 +
    1.52 +        Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
    1.53 +        COLAMD is also available under alternate licenses, contact T. Davis
    1.54 +        for details.
    1.55 +
    1.56 +        This library is free software; you can redistribute it and/or
    1.57 +        modify it under the terms of the GNU Lesser General Public
    1.58 +        License as published by the Free Software Foundation; either
    1.59 +        version 2.1 of the License, or (at your option) any later version.
    1.60 +
    1.61 +        This library is distributed in the hope that it will be useful,
    1.62 +        but WITHOUT ANY WARRANTY; without even the implied warranty of
    1.63 +        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    1.64 +        Lesser General Public License for more details.
    1.65 +
    1.66 +        You should have received a copy of the GNU Lesser General Public
    1.67 +        License along with this library; if not, write to the Free Software
    1.68 +        Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
    1.69 +        USA
    1.70 +
    1.71 +        Permission is hereby granted to use or copy this program under the
    1.72 +        terms of the GNU LGPL, provided that the Copyright, this License,
    1.73 +        and the Availability of the original version is retained on all copies.
    1.74 +        User documentation of any code that uses this code or any modified
    1.75 +        version of this code must cite the Copyright, this License, the
    1.76 +        Availability note, and "Used by permission." Permission to modify
    1.77 +        the code and to distribute modified code is granted, provided the
    1.78 +        Copyright, this License, and the Availability note are retained,
    1.79 +        and a notice that the code was modified is included.
    1.80 +
    1.81 +    Availability:
    1.82 +
    1.83 +        The colamd/symamd library is available at
    1.84 +
    1.85 +            http://www.cise.ufl.edu/research/sparse/colamd/
    1.86 +
    1.87 +        This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
    1.88 +        file.  It requires the colamd.h file.  It is required by the colamdmex.c
    1.89 +        and symamdmex.c files, for the MATLAB interface to colamd and symamd.
    1.90 +        Appears as ACM Algorithm 836.
    1.91 +
    1.92 +    See the ChangeLog file for changes since Version 1.0.
    1.93 +
    1.94 +    References:
    1.95 +
    1.96 +        T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
    1.97 +        minimum degree ordering algorithm, ACM Transactions on Mathematical
    1.98 +        Software, vol. 30, no. 3., pp. 353-376, 2004.
    1.99 +
   1.100 +        T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
   1.101 +        an approximate column minimum degree ordering algorithm, ACM
   1.102 +        Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
   1.103 +        2004.
   1.104 +
   1.105 +*/
   1.106 +
   1.107 +/* ========================================================================== */
   1.108 +/* === Description of user-callable routines ================================ */
   1.109 +/* ========================================================================== */
   1.110 +
   1.111 +/* COLAMD includes both int and UF_long versions of all its routines.  The
   1.112 + * description below is for the int version.  For UF_long, all int arguments
   1.113 + * become UF_long.  UF_long is normally defined as long, except for WIN64.
   1.114 +
   1.115 +    ----------------------------------------------------------------------------
   1.116 +    colamd_recommended:
   1.117 +    ----------------------------------------------------------------------------
   1.118 +
   1.119 +        C syntax:
   1.120 +
   1.121 +            #include "colamd.h"
   1.122 +            size_t colamd_recommended (int nnz, int n_row, int n_col) ;
   1.123 +            size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
   1.124 +                UF_long n_col) ;
   1.125 +
   1.126 +        Purpose:
   1.127 +
   1.128 +            Returns recommended value of Alen for use by colamd.  Returns 0
   1.129 +            if any input argument is negative.  The use of this routine
   1.130 +            is optional.  Not needed for symamd, which dynamically allocates
   1.131 +            its own memory.
   1.132 +
   1.133 +            Note that in v2.4 and earlier, these routines returned int or long.
   1.134 +            They now return a value of type size_t.
   1.135 +
   1.136 +        Arguments (all input arguments):
   1.137 +
   1.138 +            int nnz ;           Number of nonzeros in the matrix A.  This must
   1.139 +                                be the same value as p [n_col] in the call to
   1.140 +                                colamd - otherwise you will get a wrong value
   1.141 +                                of the recommended memory to use.
   1.142 +
   1.143 +            int n_row ;         Number of rows in the matrix A.
   1.144 +
   1.145 +            int n_col ;         Number of columns in the matrix A.
   1.146 +
   1.147 +    ----------------------------------------------------------------------------
   1.148 +    colamd_set_defaults:
   1.149 +    ----------------------------------------------------------------------------
   1.150 +
   1.151 +        C syntax:
   1.152 +
   1.153 +            #include "colamd.h"
   1.154 +            colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
   1.155 +            colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
   1.156 +
   1.157 +        Purpose:
   1.158 +
   1.159 +            Sets the default parameters.  The use of this routine is optional.
   1.160 +
   1.161 +        Arguments:
   1.162 +
   1.163 +            double knobs [COLAMD_KNOBS] ;       Output only.
   1.164 +
   1.165 +                NOTE: the meaning of the dense row/col knobs has changed in v2.4
   1.166 +
   1.167 +                knobs [0] and knobs [1] control dense row and col detection:
   1.168 +
   1.169 +                Colamd: rows with more than
   1.170 +                max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
   1.171 +                entries are removed prior to ordering.  Columns with more than
   1.172 +                max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
   1.173 +                entries are removed prior to
   1.174 +                ordering, and placed last in the output column ordering. 
   1.175 +
   1.176 +                Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
   1.177 +                Rows and columns with more than
   1.178 +                max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
   1.179 +                entries are removed prior to ordering, and placed last in the
   1.180 +                output ordering.
   1.181 +
   1.182 +                COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
   1.183 +                respectively, in colamd.h.  Default values of these two knobs
   1.184 +                are both 10.  Currently, only knobs [0] and knobs [1] are
   1.185 +                used, but future versions may use more knobs.  If so, they will
   1.186 +                be properly set to their defaults by the future version of
   1.187 +                colamd_set_defaults, so that the code that calls colamd will
   1.188 +                not need to change, assuming that you either use
   1.189 +                colamd_set_defaults, or pass a (double *) NULL pointer as the
   1.190 +                knobs array to colamd or symamd.
   1.191 +
   1.192 +            knobs [2]: aggressive absorption
   1.193 +
   1.194 +                knobs [COLAMD_AGGRESSIVE] controls whether or not to do
   1.195 +                aggressive absorption during the ordering.  Default is TRUE.
   1.196 +
   1.197 +
   1.198 +    ----------------------------------------------------------------------------
   1.199 +    colamd:
   1.200 +    ----------------------------------------------------------------------------
   1.201 +
   1.202 +        C syntax:
   1.203 +
   1.204 +            #include "colamd.h"
   1.205 +            int colamd (int n_row, int n_col, int Alen, int *A, int *p,
   1.206 +                double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
   1.207 +            UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
   1.208 +                UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
   1.209 +                UF_long stats [COLAMD_STATS]) ;
   1.210 +
   1.211 +        Purpose:
   1.212 +
   1.213 +            Computes a column ordering (Q) of A such that P(AQ)=LU or
   1.214 +            (AQ)'AQ=LL' have less fill-in and require fewer floating point
   1.215 +            operations than factorizing the unpermuted matrix A or A'A,
   1.216 +            respectively.
   1.217 +            
   1.218 +        Returns:
   1.219 +
   1.220 +            TRUE (1) if successful, FALSE (0) otherwise.
   1.221 +
   1.222 +        Arguments:
   1.223 +
   1.224 +            int n_row ;         Input argument.
   1.225 +
   1.226 +                Number of rows in the matrix A.
   1.227 +                Restriction:  n_row >= 0.
   1.228 +                Colamd returns FALSE if n_row is negative.
   1.229 +
   1.230 +            int n_col ;         Input argument.
   1.231 +
   1.232 +                Number of columns in the matrix A.
   1.233 +                Restriction:  n_col >= 0.
   1.234 +                Colamd returns FALSE if n_col is negative.
   1.235 +
   1.236 +            int Alen ;          Input argument.
   1.237 +
   1.238 +                Restriction (see note):
   1.239 +                Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
   1.240 +                Colamd returns FALSE if these conditions are not met.
   1.241 +
   1.242 +                Note:  this restriction makes an modest assumption regarding
   1.243 +                the size of the two typedef's structures in colamd.h.
   1.244 +                We do, however, guarantee that
   1.245 +
   1.246 +                        Alen >= colamd_recommended (nnz, n_row, n_col)
   1.247 +
   1.248 +                will be sufficient.  Note: the macro version does not check
   1.249 +                for integer overflow, and thus is not recommended.  Use
   1.250 +                the colamd_recommended routine instead.
   1.251 +
   1.252 +            int A [Alen] ;      Input argument, undefined on output.
   1.253 +
   1.254 +                A is an integer array of size Alen.  Alen must be at least as
   1.255 +                large as the bare minimum value given above, but this is very
   1.256 +                low, and can result in excessive run time.  For best
   1.257 +                performance, we recommend that Alen be greater than or equal to
   1.258 +                colamd_recommended (nnz, n_row, n_col), which adds
   1.259 +                nnz/5 to the bare minimum value given above.
   1.260 +
   1.261 +                On input, the row indices of the entries in column c of the
   1.262 +                matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
   1.263 +                in a given column c need not be in ascending order, and
   1.264 +                duplicate row indices may be be present.  However, colamd will
   1.265 +                work a little faster if both of these conditions are met
   1.266 +                (Colamd puts the matrix into this format, if it finds that the
   1.267 +                the conditions are not met).
   1.268 +
   1.269 +                The matrix is 0-based.  That is, rows are in the range 0 to
   1.270 +                n_row-1, and columns are in the range 0 to n_col-1.  Colamd
   1.271 +                returns FALSE if any row index is out of range.
   1.272 +
   1.273 +                The contents of A are modified during ordering, and are
   1.274 +                undefined on output.
   1.275 +
   1.276 +            int p [n_col+1] ;   Both input and output argument.
   1.277 +
   1.278 +                p is an integer array of size n_col+1.  On input, it holds the
   1.279 +                "pointers" for the column form of the matrix A.  Column c of
   1.280 +                the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
   1.281 +                entry, p [0], must be zero, and p [c] <= p [c+1] must hold
   1.282 +                for all c in the range 0 to n_col-1.  The value p [n_col] is
   1.283 +                thus the total number of entries in the pattern of the matrix A.
   1.284 +                Colamd returns FALSE if these conditions are not met.
   1.285 +
   1.286 +                On output, if colamd returns TRUE, the array p holds the column
   1.287 +                permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
   1.288 +                the first column index in the new ordering, and p [n_col-1] is
   1.289 +                the last.  That is, p [k] = j means that column j of A is the
   1.290 +                kth pivot column, in AQ, where k is in the range 0 to n_col-1
   1.291 +                (p [0] = j means that column j of A is the first column in AQ).
   1.292 +
   1.293 +                If colamd returns FALSE, then no permutation is returned, and
   1.294 +                p is undefined on output.
   1.295 +
   1.296 +            double knobs [COLAMD_KNOBS] ;       Input argument.
   1.297 +
   1.298 +                See colamd_set_defaults for a description.
   1.299 +
   1.300 +            int stats [COLAMD_STATS] ;          Output argument.
   1.301 +
   1.302 +                Statistics on the ordering, and error status.
   1.303 +                See colamd.h for related definitions.
   1.304 +                Colamd returns FALSE if stats is not present.
   1.305 +
   1.306 +                stats [0]:  number of dense or empty rows ignored.
   1.307 +
   1.308 +                stats [1]:  number of dense or empty columns ignored (and
   1.309 +                                ordered last in the output permutation p)
   1.310 +                                Note that a row can become "empty" if it
   1.311 +                                contains only "dense" and/or "empty" columns,
   1.312 +                                and similarly a column can become "empty" if it
   1.313 +                                only contains "dense" and/or "empty" rows.
   1.314 +
   1.315 +                stats [2]:  number of garbage collections performed.
   1.316 +                                This can be excessively high if Alen is close
   1.317 +                                to the minimum required value.
   1.318 +
   1.319 +                stats [3]:  status code.  < 0 is an error code.
   1.320 +                            > 1 is a warning or notice.
   1.321 +
   1.322 +                        0       OK.  Each column of the input matrix contained
   1.323 +                                row indices in increasing order, with no
   1.324 +                                duplicates.
   1.325 +
   1.326 +                        1       OK, but columns of input matrix were jumbled
   1.327 +                                (unsorted columns or duplicate entries).  Colamd
   1.328 +                                had to do some extra work to sort the matrix
   1.329 +                                first and remove duplicate entries, but it
   1.330 +                                still was able to return a valid permutation
   1.331 +                                (return value of colamd was TRUE).
   1.332 +
   1.333 +                                        stats [4]: highest numbered column that
   1.334 +                                                is unsorted or has duplicate
   1.335 +                                                entries.
   1.336 +                                        stats [5]: last seen duplicate or
   1.337 +                                                unsorted row index.
   1.338 +                                        stats [6]: number of duplicate or
   1.339 +                                                unsorted row indices.
   1.340 +
   1.341 +                        -1      A is a null pointer
   1.342 +
   1.343 +                        -2      p is a null pointer
   1.344 +
   1.345 +                        -3      n_row is negative
   1.346 +
   1.347 +                                        stats [4]: n_row
   1.348 +
   1.349 +                        -4      n_col is negative
   1.350 +
   1.351 +                                        stats [4]: n_col
   1.352 +
   1.353 +                        -5      number of nonzeros in matrix is negative
   1.354 +
   1.355 +                                        stats [4]: number of nonzeros, p [n_col]
   1.356 +
   1.357 +                        -6      p [0] is nonzero
   1.358 +
   1.359 +                                        stats [4]: p [0]
   1.360 +
   1.361 +                        -7      A is too small
   1.362 +
   1.363 +                                        stats [4]: required size
   1.364 +                                        stats [5]: actual size (Alen)
   1.365 +
   1.366 +                        -8      a column has a negative number of entries
   1.367 +
   1.368 +                                        stats [4]: column with < 0 entries
   1.369 +                                        stats [5]: number of entries in col
   1.370 +
   1.371 +                        -9      a row index is out of bounds
   1.372 +
   1.373 +                                        stats [4]: column with bad row index
   1.374 +                                        stats [5]: bad row index
   1.375 +                                        stats [6]: n_row, # of rows of matrx
   1.376 +
   1.377 +                        -10     (unused; see symamd.c)
   1.378 +
   1.379 +                        -999    (unused; see symamd.c)
   1.380 +
   1.381 +                Future versions may return more statistics in the stats array.
   1.382 +
   1.383 +        Example:
   1.384 +        
   1.385 +            See http://www.cise.ufl.edu/research/sparse/colamd/example.c
   1.386 +            for a complete example.
   1.387 +
   1.388 +            To order the columns of a 5-by-4 matrix with 11 nonzero entries in
   1.389 +            the following nonzero pattern
   1.390 +
   1.391 +                x 0 x 0
   1.392 +                x 0 x x
   1.393 +                0 x x 0
   1.394 +                0 0 x x
   1.395 +                x x 0 0
   1.396 +
   1.397 +            with default knobs and no output statistics, do the following:
   1.398 +
   1.399 +                #include "colamd.h"
   1.400 +                #define ALEN 100
   1.401 +                int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
   1.402 +                int p [ ] = {0, 3, 5, 9, 11} ;
   1.403 +                int stats [COLAMD_STATS] ;
   1.404 +                colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
   1.405 +
   1.406 +            The permutation is returned in the array p, and A is destroyed.
   1.407 +
   1.408 +    ----------------------------------------------------------------------------
   1.409 +    symamd:
   1.410 +    ----------------------------------------------------------------------------
   1.411 +
   1.412 +        C syntax:
   1.413 +
   1.414 +            #include "colamd.h"
   1.415 +            int symamd (int n, int *A, int *p, int *perm,
   1.416 +                double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
   1.417 +                void (*allocate) (size_t, size_t), void (*release) (void *)) ;
   1.418 +            UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
   1.419 +                double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
   1.420 +                void (*allocate) (size_t, size_t), void (*release) (void *)) ;
   1.421 +
   1.422 +        Purpose:
   1.423 +
   1.424 +            The symamd routine computes an ordering P of a symmetric sparse
   1.425 +            matrix A such that the Cholesky factorization PAP' = LL' remains
   1.426 +            sparse.  It is based on a column ordering of a matrix M constructed
   1.427 +            so that the nonzero pattern of M'M is the same as A.  The matrix A
   1.428 +            is assumed to be symmetric; only the strictly lower triangular part
   1.429 +            is accessed.  You must pass your selected memory allocator (usually
   1.430 +            calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
   1.431 +            memory for the temporary matrix M.
   1.432 +
   1.433 +        Returns:
   1.434 +
   1.435 +            TRUE (1) if successful, FALSE (0) otherwise.
   1.436 +
   1.437 +        Arguments:
   1.438 +
   1.439 +            int n ;             Input argument.
   1.440 +
   1.441 +                Number of rows and columns in the symmetrix matrix A.
   1.442 +                Restriction:  n >= 0.
   1.443 +                Symamd returns FALSE if n is negative.
   1.444 +
   1.445 +            int A [nnz] ;       Input argument.
   1.446 +
   1.447 +                A is an integer array of size nnz, where nnz = p [n].
   1.448 +                
   1.449 +                The row indices of the entries in column c of the matrix are
   1.450 +                held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
   1.451 +                given column c need not be in ascending order, and duplicate
   1.452 +                row indices may be present.  However, symamd will run faster
   1.453 +                if the columns are in sorted order with no duplicate entries. 
   1.454 +
   1.455 +                The matrix is 0-based.  That is, rows are in the range 0 to
   1.456 +                n-1, and columns are in the range 0 to n-1.  Symamd
   1.457 +                returns FALSE if any row index is out of range.
   1.458 +
   1.459 +                The contents of A are not modified.
   1.460 +
   1.461 +            int p [n+1] ;       Input argument.
   1.462 +
   1.463 +                p is an integer array of size n+1.  On input, it holds the
   1.464 +                "pointers" for the column form of the matrix A.  Column c of
   1.465 +                the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
   1.466 +                entry, p [0], must be zero, and p [c] <= p [c+1] must hold
   1.467 +                for all c in the range 0 to n-1.  The value p [n] is
   1.468 +                thus the total number of entries in the pattern of the matrix A.
   1.469 +                Symamd returns FALSE if these conditions are not met.
   1.470 +
   1.471 +                The contents of p are not modified.
   1.472 +
   1.473 +            int perm [n+1] ;    Output argument.
   1.474 +
   1.475 +                On output, if symamd returns TRUE, the array perm holds the
   1.476 +                permutation P, where perm [0] is the first index in the new
   1.477 +                ordering, and perm [n-1] is the last.  That is, perm [k] = j
   1.478 +                means that row and column j of A is the kth column in PAP',
   1.479 +                where k is in the range 0 to n-1 (perm [0] = j means
   1.480 +                that row and column j of A are the first row and column in
   1.481 +                PAP').  The array is used as a workspace during the ordering,
   1.482 +                which is why it must be of length n+1, not just n.
   1.483 +
   1.484 +            double knobs [COLAMD_KNOBS] ;       Input argument.
   1.485 +
   1.486 +                See colamd_set_defaults for a description.
   1.487 +
   1.488 +            int stats [COLAMD_STATS] ;          Output argument.
   1.489 +
   1.490 +                Statistics on the ordering, and error status.
   1.491 +                See colamd.h for related definitions.
   1.492 +                Symamd returns FALSE if stats is not present.
   1.493 +
   1.494 +                stats [0]:  number of dense or empty row and columns ignored
   1.495 +                                (and ordered last in the output permutation 
   1.496 +                                perm).  Note that a row/column can become
   1.497 +                                "empty" if it contains only "dense" and/or
   1.498 +                                "empty" columns/rows.
   1.499 +
   1.500 +                stats [1]:  (same as stats [0])
   1.501 +
   1.502 +                stats [2]:  number of garbage collections performed.
   1.503 +
   1.504 +                stats [3]:  status code.  < 0 is an error code.
   1.505 +                            > 1 is a warning or notice.
   1.506 +
   1.507 +                        0       OK.  Each column of the input matrix contained
   1.508 +                                row indices in increasing order, with no
   1.509 +                                duplicates.
   1.510 +
   1.511 +                        1       OK, but columns of input matrix were jumbled
   1.512 +                                (unsorted columns or duplicate entries).  Symamd
   1.513 +                                had to do some extra work to sort the matrix
   1.514 +                                first and remove duplicate entries, but it
   1.515 +                                still was able to return a valid permutation
   1.516 +                                (return value of symamd was TRUE).
   1.517 +
   1.518 +                                        stats [4]: highest numbered column that
   1.519 +                                                is unsorted or has duplicate
   1.520 +                                                entries.
   1.521 +                                        stats [5]: last seen duplicate or
   1.522 +                                                unsorted row index.
   1.523 +                                        stats [6]: number of duplicate or
   1.524 +                                                unsorted row indices.
   1.525 +
   1.526 +                        -1      A is a null pointer
   1.527 +
   1.528 +                        -2      p is a null pointer
   1.529 +
   1.530 +                        -3      (unused, see colamd.c)
   1.531 +
   1.532 +                        -4      n is negative
   1.533 +
   1.534 +                                        stats [4]: n
   1.535 +
   1.536 +                        -5      number of nonzeros in matrix is negative
   1.537 +
   1.538 +                                        stats [4]: # of nonzeros (p [n]).
   1.539 +
   1.540 +                        -6      p [0] is nonzero
   1.541 +
   1.542 +                                        stats [4]: p [0]
   1.543 +
   1.544 +                        -7      (unused)
   1.545 +
   1.546 +                        -8      a column has a negative number of entries
   1.547 +
   1.548 +                                        stats [4]: column with < 0 entries
   1.549 +                                        stats [5]: number of entries in col
   1.550 +
   1.551 +                        -9      a row index is out of bounds
   1.552 +
   1.553 +                                        stats [4]: column with bad row index
   1.554 +                                        stats [5]: bad row index
   1.555 +                                        stats [6]: n_row, # of rows of matrx
   1.556 +
   1.557 +                        -10     out of memory (unable to allocate temporary
   1.558 +                                workspace for M or count arrays using the
   1.559 +                                "allocate" routine passed into symamd).
   1.560 +
   1.561 +                Future versions may return more statistics in the stats array.
   1.562 +
   1.563 +            void * (*allocate) (size_t, size_t)
   1.564 +
   1.565 +                A pointer to a function providing memory allocation.  The
   1.566 +                allocated memory must be returned initialized to zero.  For a
   1.567 +                C application, this argument should normally be a pointer to
   1.568 +                calloc.  For a MATLAB mexFunction, the routine mxCalloc is
   1.569 +                passed instead.
   1.570 +
   1.571 +            void (*release) (size_t, size_t)
   1.572 +
   1.573 +                A pointer to a function that frees memory allocated by the
   1.574 +                memory allocation routine above.  For a C application, this
   1.575 +                argument should normally be a pointer to free.  For a MATLAB
   1.576 +                mexFunction, the routine mxFree is passed instead.
   1.577 +
   1.578 +
   1.579 +    ----------------------------------------------------------------------------
   1.580 +    colamd_report:
   1.581 +    ----------------------------------------------------------------------------
   1.582 +
   1.583 +        C syntax:
   1.584 +
   1.585 +            #include "colamd.h"
   1.586 +            colamd_report (int stats [COLAMD_STATS]) ;
   1.587 +            colamd_l_report (UF_long stats [COLAMD_STATS]) ;
   1.588 +
   1.589 +        Purpose:
   1.590 +
   1.591 +            Prints the error status and statistics recorded in the stats
   1.592 +            array on the standard error output (for a standard C routine)
   1.593 +            or on the MATLAB output (for a mexFunction).
   1.594 +
   1.595 +        Arguments:
   1.596 +
   1.597 +            int stats [COLAMD_STATS] ;  Input only.  Statistics from colamd.
   1.598 +
   1.599 +
   1.600 +    ----------------------------------------------------------------------------
   1.601 +    symamd_report:
   1.602 +    ----------------------------------------------------------------------------
   1.603 +
   1.604 +        C syntax:
   1.605 +
   1.606 +            #include "colamd.h"
   1.607 +            symamd_report (int stats [COLAMD_STATS]) ;
   1.608 +            symamd_l_report (UF_long stats [COLAMD_STATS]) ;
   1.609 +
   1.610 +        Purpose:
   1.611 +
   1.612 +            Prints the error status and statistics recorded in the stats
   1.613 +            array on the standard error output (for a standard C routine)
   1.614 +            or on the MATLAB output (for a mexFunction).
   1.615 +
   1.616 +        Arguments:
   1.617 +
   1.618 +            int stats [COLAMD_STATS] ;  Input only.  Statistics from symamd.
   1.619 +
   1.620 +
   1.621 +*/
   1.622 +
   1.623 +/* ========================================================================== */
   1.624 +/* === Scaffolding code definitions  ======================================== */
   1.625 +/* ========================================================================== */
   1.626 +
   1.627 +/* Ensure that debugging is turned off: */
   1.628 +#ifndef NDEBUG
   1.629 +#define NDEBUG
   1.630 +#endif
   1.631 +
   1.632 +/* turn on debugging by uncommenting the following line
   1.633 + #undef NDEBUG
   1.634 +*/
   1.635 +
   1.636 +/*
   1.637 +   Our "scaffolding code" philosophy:  In our opinion, well-written library
   1.638 +   code should keep its "debugging" code, and just normally have it turned off
   1.639 +   by the compiler so as not to interfere with performance.  This serves
   1.640 +   several purposes:
   1.641 +
   1.642 +   (1) assertions act as comments to the reader, telling you what the code
   1.643 +        expects at that point.  All assertions will always be true (unless
   1.644 +        there really is a bug, of course).
   1.645 +
   1.646 +   (2) leaving in the scaffolding code assists anyone who would like to modify
   1.647 +        the code, or understand the algorithm (by reading the debugging output,
   1.648 +        one can get a glimpse into what the code is doing).
   1.649 +
   1.650 +   (3) (gasp!) for actually finding bugs.  This code has been heavily tested
   1.651 +        and "should" be fully functional and bug-free ... but you never know...
   1.652 +
   1.653 +    The code will become outrageously slow when debugging is
   1.654 +    enabled.  To control the level of debugging output, set an environment
   1.655 +    variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
   1.656 +    you should see the following message on the standard output:
   1.657 +
   1.658 +        colamd: debug version, D = 1 (THIS WILL BE SLOW!)
   1.659 +
   1.660 +    or a similar message for symamd.  If you don't, then debugging has not
   1.661 +    been enabled.
   1.662 +
   1.663 +*/
   1.664 +
   1.665 +/* ========================================================================== */
   1.666 +/* === Include files ======================================================== */
   1.667 +/* ========================================================================== */
   1.668 +
   1.669 +#include "colamd.h"
   1.670 +
   1.671 +#if 0 /* by mao */
   1.672 +#include <limits.h>
   1.673 +#include <math.h>
   1.674 +
   1.675 +#ifdef MATLAB_MEX_FILE
   1.676 +#include "mex.h"
   1.677 +#include "matrix.h"
   1.678 +#endif /* MATLAB_MEX_FILE */
   1.679 +
   1.680 +#if !defined (NPRINT) || !defined (NDEBUG)
   1.681 +#include <stdio.h>
   1.682 +#endif
   1.683 +
   1.684 +#ifndef NULL
   1.685 +#define NULL ((void *) 0)
   1.686 +#endif
   1.687 +#endif
   1.688 +
   1.689 +/* ========================================================================== */
   1.690 +/* === int or UF_long ======================================================= */
   1.691 +/* ========================================================================== */
   1.692 +
   1.693 +#if 0 /* by mao */
   1.694 +/* define UF_long */
   1.695 +#include "UFconfig.h"
   1.696 +#endif
   1.697 +
   1.698 +#ifdef DLONG
   1.699 +
   1.700 +#define Int UF_long
   1.701 +#define ID  UF_long_id
   1.702 +#define Int_MAX UF_long_max
   1.703 +
   1.704 +#define COLAMD_recommended colamd_l_recommended
   1.705 +#define COLAMD_set_defaults colamd_l_set_defaults
   1.706 +#define COLAMD_MAIN colamd_l
   1.707 +#define SYMAMD_MAIN symamd_l
   1.708 +#define COLAMD_report colamd_l_report
   1.709 +#define SYMAMD_report symamd_l_report
   1.710 +
   1.711 +#else
   1.712 +
   1.713 +#define Int int
   1.714 +#define ID "%d"
   1.715 +#define Int_MAX INT_MAX
   1.716 +
   1.717 +#define COLAMD_recommended colamd_recommended
   1.718 +#define COLAMD_set_defaults colamd_set_defaults
   1.719 +#define COLAMD_MAIN colamd
   1.720 +#define SYMAMD_MAIN symamd
   1.721 +#define COLAMD_report colamd_report
   1.722 +#define SYMAMD_report symamd_report
   1.723 +
   1.724 +#endif
   1.725 +
   1.726 +/* ========================================================================== */
   1.727 +/* === Row and Column structures ============================================ */
   1.728 +/* ========================================================================== */
   1.729 +
   1.730 +/* User code that makes use of the colamd/symamd routines need not directly */
   1.731 +/* reference these structures.  They are used only for colamd_recommended. */
   1.732 +
   1.733 +typedef struct Colamd_Col_struct
   1.734 +{
   1.735 +    Int start ;         /* index for A of first row in this column, or DEAD */
   1.736 +                        /* if column is dead */
   1.737 +    Int length ;        /* number of rows in this column */
   1.738 +    union
   1.739 +    {
   1.740 +        Int thickness ; /* number of original columns represented by this */
   1.741 +                        /* col, if the column is alive */
   1.742 +        Int parent ;    /* parent in parent tree super-column structure, if */
   1.743 +                        /* the column is dead */
   1.744 +    } shared1 ;
   1.745 +    union
   1.746 +    {
   1.747 +        Int score ;     /* the score used to maintain heap, if col is alive */
   1.748 +        Int order ;     /* pivot ordering of this column, if col is dead */
   1.749 +    } shared2 ;
   1.750 +    union
   1.751 +    {
   1.752 +        Int headhash ;  /* head of a hash bucket, if col is at the head of */
   1.753 +                        /* a degree list */
   1.754 +        Int hash ;      /* hash value, if col is not in a degree list */
   1.755 +        Int prev ;      /* previous column in degree list, if col is in a */
   1.756 +                        /* degree list (but not at the head of a degree list) */
   1.757 +    } shared3 ;
   1.758 +    union
   1.759 +    {
   1.760 +        Int degree_next ;       /* next column, if col is in a degree list */
   1.761 +        Int hash_next ;         /* next column, if col is in a hash list */
   1.762 +    } shared4 ;
   1.763 +
   1.764 +} Colamd_Col ;
   1.765 +
   1.766 +typedef struct Colamd_Row_struct
   1.767 +{
   1.768 +    Int start ;         /* index for A of first col in this row */
   1.769 +    Int length ;        /* number of principal columns in this row */
   1.770 +    union
   1.771 +    {
   1.772 +        Int degree ;    /* number of principal & non-principal columns in row */
   1.773 +        Int p ;         /* used as a row pointer in init_rows_cols () */
   1.774 +    } shared1 ;
   1.775 +    union
   1.776 +    {
   1.777 +        Int mark ;      /* for computing set differences and marking dead rows*/
   1.778 +        Int first_column ;/* first column in row (used in garbage collection) */
   1.779 +    } shared2 ;
   1.780 +
   1.781 +} Colamd_Row ;
   1.782 +
   1.783 +/* ========================================================================== */
   1.784 +/* === Definitions ========================================================== */
   1.785 +/* ========================================================================== */
   1.786 +
   1.787 +/* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
   1.788 +#define PUBLIC
   1.789 +#define PRIVATE static
   1.790 +
   1.791 +#define DENSE_DEGREE(alpha,n) \
   1.792 +    ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
   1.793 +
   1.794 +#define MAX(a,b) (((a) > (b)) ? (a) : (b))
   1.795 +#define MIN(a,b) (((a) < (b)) ? (a) : (b))
   1.796 +
   1.797 +#define ONES_COMPLEMENT(r) (-(r)-1)
   1.798 +
   1.799 +/* -------------------------------------------------------------------------- */
   1.800 +/* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
   1.801 +/* -------------------------------------------------------------------------- */
   1.802 +
   1.803 +#ifndef TRUE
   1.804 +#define TRUE (1)
   1.805 +#endif
   1.806 +
   1.807 +#ifndef FALSE
   1.808 +#define FALSE (0)
   1.809 +#endif
   1.810 +
   1.811 +/* -------------------------------------------------------------------------- */
   1.812 +
   1.813 +#define EMPTY   (-1)
   1.814 +
   1.815 +/* Row and column status */
   1.816 +#define ALIVE   (0)
   1.817 +#define DEAD    (-1)
   1.818 +
   1.819 +/* Column status */
   1.820 +#define DEAD_PRINCIPAL          (-1)
   1.821 +#define DEAD_NON_PRINCIPAL      (-2)
   1.822 +
   1.823 +/* Macros for row and column status update and checking. */
   1.824 +#define ROW_IS_DEAD(r)                  ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
   1.825 +#define ROW_IS_MARKED_DEAD(row_mark)    (row_mark < ALIVE)
   1.826 +#define ROW_IS_ALIVE(r)                 (Row [r].shared2.mark >= ALIVE)
   1.827 +#define COL_IS_DEAD(c)                  (Col [c].start < ALIVE)
   1.828 +#define COL_IS_ALIVE(c)                 (Col [c].start >= ALIVE)
   1.829 +#define COL_IS_DEAD_PRINCIPAL(c)        (Col [c].start == DEAD_PRINCIPAL)
   1.830 +#define KILL_ROW(r)                     { Row [r].shared2.mark = DEAD ; }
   1.831 +#define KILL_PRINCIPAL_COL(c)           { Col [c].start = DEAD_PRINCIPAL ; }
   1.832 +#define KILL_NON_PRINCIPAL_COL(c)       { Col [c].start = DEAD_NON_PRINCIPAL ; }
   1.833 +
   1.834 +/* ========================================================================== */
   1.835 +/* === Colamd reporting mechanism =========================================== */
   1.836 +/* ========================================================================== */
   1.837 +
   1.838 +#if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
   1.839 +/* In MATLAB, matrices are 1-based to the user, but 0-based internally */
   1.840 +#define INDEX(i) ((i)+1)
   1.841 +#else
   1.842 +/* In C, matrices are 0-based and indices are reported as such in *_report */
   1.843 +#define INDEX(i) (i)
   1.844 +#endif
   1.845 +
   1.846 +/* All output goes through the PRINTF macro.  */
   1.847 +#define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
   1.848 +
   1.849 +/* ========================================================================== */
   1.850 +/* === Prototypes of PRIVATE routines ======================================= */
   1.851 +/* ========================================================================== */
   1.852 +
   1.853 +PRIVATE Int init_rows_cols
   1.854 +(
   1.855 +    Int n_row,
   1.856 +    Int n_col,
   1.857 +    Colamd_Row Row [],
   1.858 +    Colamd_Col Col [],
   1.859 +    Int A [],
   1.860 +    Int p [],
   1.861 +    Int stats [COLAMD_STATS]
   1.862 +) ;
   1.863 +
   1.864 +PRIVATE void init_scoring
   1.865 +(
   1.866 +    Int n_row,
   1.867 +    Int n_col,
   1.868 +    Colamd_Row Row [],
   1.869 +    Colamd_Col Col [],
   1.870 +    Int A [],
   1.871 +    Int head [],
   1.872 +    double knobs [COLAMD_KNOBS],
   1.873 +    Int *p_n_row2,
   1.874 +    Int *p_n_col2,
   1.875 +    Int *p_max_deg
   1.876 +) ;
   1.877 +
   1.878 +PRIVATE Int find_ordering
   1.879 +(
   1.880 +    Int n_row,
   1.881 +    Int n_col,
   1.882 +    Int Alen,
   1.883 +    Colamd_Row Row [],
   1.884 +    Colamd_Col Col [],
   1.885 +    Int A [],
   1.886 +    Int head [],
   1.887 +    Int n_col2,
   1.888 +    Int max_deg,
   1.889 +    Int pfree,
   1.890 +    Int aggressive
   1.891 +) ;
   1.892 +
   1.893 +PRIVATE void order_children
   1.894 +(
   1.895 +    Int n_col,
   1.896 +    Colamd_Col Col [],
   1.897 +    Int p []
   1.898 +) ;
   1.899 +
   1.900 +PRIVATE void detect_super_cols
   1.901 +(
   1.902 +
   1.903 +#ifndef NDEBUG
   1.904 +    Int n_col,
   1.905 +    Colamd_Row Row [],
   1.906 +#endif /* NDEBUG */
   1.907 +
   1.908 +    Colamd_Col Col [],
   1.909 +    Int A [],
   1.910 +    Int head [],
   1.911 +    Int row_start,
   1.912 +    Int row_length
   1.913 +) ;
   1.914 +
   1.915 +PRIVATE Int garbage_collection
   1.916 +(
   1.917 +    Int n_row,
   1.918 +    Int n_col,
   1.919 +    Colamd_Row Row [],
   1.920 +    Colamd_Col Col [],
   1.921 +    Int A [],
   1.922 +    Int *pfree
   1.923 +) ;
   1.924 +
   1.925 +PRIVATE Int clear_mark
   1.926 +(
   1.927 +    Int tag_mark,
   1.928 +    Int max_mark,
   1.929 +    Int n_row,
   1.930 +    Colamd_Row Row []
   1.931 +) ;
   1.932 +
   1.933 +PRIVATE void print_report
   1.934 +(
   1.935 +    char *method,
   1.936 +    Int stats [COLAMD_STATS]
   1.937 +) ;
   1.938 +
   1.939 +/* ========================================================================== */
   1.940 +/* === Debugging prototypes and definitions ================================= */
   1.941 +/* ========================================================================== */
   1.942 +
   1.943 +#ifndef NDEBUG
   1.944 +
   1.945 +#if 0 /* by mao */
   1.946 +#include <assert.h>
   1.947 +#endif
   1.948 +
   1.949 +/* colamd_debug is the *ONLY* global variable, and is only */
   1.950 +/* present when debugging */
   1.951 +
   1.952 +PRIVATE Int colamd_debug = 0 ;  /* debug print level */
   1.953 +
   1.954 +#define DEBUG0(params) { PRINTF (params) ; }
   1.955 +#define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
   1.956 +#define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
   1.957 +#define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
   1.958 +#define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
   1.959 +
   1.960 +#if 0 /* by mao */
   1.961 +#ifdef MATLAB_MEX_FILE
   1.962 +#define ASSERT(expression) (mxAssert ((expression), ""))
   1.963 +#else
   1.964 +#define ASSERT(expression) (assert (expression))
   1.965 +#endif /* MATLAB_MEX_FILE */
   1.966 +#else
   1.967 +#define ASSERT xassert
   1.968 +#endif
   1.969 +
   1.970 +PRIVATE void colamd_get_debug   /* gets the debug print level from getenv */
   1.971 +(
   1.972 +    char *method
   1.973 +) ;
   1.974 +
   1.975 +PRIVATE void debug_deg_lists
   1.976 +(
   1.977 +    Int n_row,
   1.978 +    Int n_col,
   1.979 +    Colamd_Row Row [],
   1.980 +    Colamd_Col Col [],
   1.981 +    Int head [],
   1.982 +    Int min_score,
   1.983 +    Int should,
   1.984 +    Int max_deg
   1.985 +) ;
   1.986 +
   1.987 +PRIVATE void debug_mark
   1.988 +(
   1.989 +    Int n_row,
   1.990 +    Colamd_Row Row [],
   1.991 +    Int tag_mark,
   1.992 +    Int max_mark
   1.993 +) ;
   1.994 +
   1.995 +PRIVATE void debug_matrix
   1.996 +(
   1.997 +    Int n_row,
   1.998 +    Int n_col,
   1.999 +    Colamd_Row Row [],
  1.1000 +    Colamd_Col Col [],
  1.1001 +    Int A []
  1.1002 +) ;
  1.1003 +
  1.1004 +PRIVATE void debug_structures
  1.1005 +(
  1.1006 +    Int n_row,
  1.1007 +    Int n_col,
  1.1008 +    Colamd_Row Row [],
  1.1009 +    Colamd_Col Col [],
  1.1010 +    Int A [],
  1.1011 +    Int n_col2
  1.1012 +) ;
  1.1013 +
  1.1014 +#else /* NDEBUG */
  1.1015 +
  1.1016 +/* === No debugging ========================================================= */
  1.1017 +
  1.1018 +#define DEBUG0(params) ;
  1.1019 +#define DEBUG1(params) ;
  1.1020 +#define DEBUG2(params) ;
  1.1021 +#define DEBUG3(params) ;
  1.1022 +#define DEBUG4(params) ;
  1.1023 +
  1.1024 +#define ASSERT(expression)
  1.1025 +
  1.1026 +#endif /* NDEBUG */
  1.1027 +
  1.1028 +/* ========================================================================== */
  1.1029 +/* === USER-CALLABLE ROUTINES: ============================================== */
  1.1030 +/* ========================================================================== */
  1.1031 +
  1.1032 +/* ========================================================================== */
  1.1033 +/* === colamd_recommended =================================================== */
  1.1034 +/* ========================================================================== */
  1.1035 +
  1.1036 +/*
  1.1037 +    The colamd_recommended routine returns the suggested size for Alen.  This
  1.1038 +    value has been determined to provide good balance between the number of
  1.1039 +    garbage collections and the memory requirements for colamd.  If any
  1.1040 +    argument is negative, or if integer overflow occurs, a 0 is returned as an
  1.1041 +    error condition.  2*nnz space is required for the row and column
  1.1042 +    indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
  1.1043 +    required for the Col and Row arrays, respectively, which are internal to
  1.1044 +    colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the
  1.1045 +    minimal amount of "elbow room", and nnz/5 more space is recommended for
  1.1046 +    run time efficiency.
  1.1047 +
  1.1048 +    Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
  1.1049 +
  1.1050 +    This function is not needed when using symamd.
  1.1051 +*/
  1.1052 +
  1.1053 +/* add two values of type size_t, and check for integer overflow */
  1.1054 +static size_t t_add (size_t a, size_t b, int *ok)
  1.1055 +{
  1.1056 +    (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
  1.1057 +    return ((*ok) ? (a + b) : 0) ;
  1.1058 +}
  1.1059 +
  1.1060 +/* compute a*k where k is a small integer, and check for integer overflow */
  1.1061 +static size_t t_mult (size_t a, size_t k, int *ok)
  1.1062 +{
  1.1063 +    size_t i, s = 0 ;
  1.1064 +    for (i = 0 ; i < k ; i++)
  1.1065 +    {
  1.1066 +        s = t_add (s, a, ok) ;
  1.1067 +    }
  1.1068 +    return (s) ;
  1.1069 +}
  1.1070 +
  1.1071 +/* size of the Col and Row structures */
  1.1072 +#define COLAMD_C(n_col,ok) \
  1.1073 +    ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
  1.1074 +
  1.1075 +#define COLAMD_R(n_row,ok) \
  1.1076 +    ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
  1.1077 +
  1.1078 +
  1.1079 +PUBLIC size_t COLAMD_recommended        /* returns recommended value of Alen. */
  1.1080 +(
  1.1081 +    /* === Parameters ======================================================= */
  1.1082 +
  1.1083 +    Int nnz,                    /* number of nonzeros in A */
  1.1084 +    Int n_row,                  /* number of rows in A */
  1.1085 +    Int n_col                   /* number of columns in A */
  1.1086 +)
  1.1087 +{
  1.1088 +    size_t s, c, r ;
  1.1089 +    int ok = TRUE ;
  1.1090 +    if (nnz < 0 || n_row < 0 || n_col < 0)
  1.1091 +    {
  1.1092 +        return (0) ;
  1.1093 +    }
  1.1094 +    s = t_mult (nnz, 2, &ok) ;      /* 2*nnz */
  1.1095 +    c = COLAMD_C (n_col, &ok) ;     /* size of column structures */
  1.1096 +    r = COLAMD_R (n_row, &ok) ;     /* size of row structures */
  1.1097 +    s = t_add (s, c, &ok) ;
  1.1098 +    s = t_add (s, r, &ok) ;
  1.1099 +    s = t_add (s, n_col, &ok) ;     /* elbow room */
  1.1100 +    s = t_add (s, nnz/5, &ok) ;     /* elbow room */
  1.1101 +    ok = ok && (s < Int_MAX) ;
  1.1102 +    return (ok ? s : 0) ;
  1.1103 +}
  1.1104 +
  1.1105 +
  1.1106 +/* ========================================================================== */
  1.1107 +/* === colamd_set_defaults ================================================== */
  1.1108 +/* ========================================================================== */
  1.1109 +
  1.1110 +/*
  1.1111 +    The colamd_set_defaults routine sets the default values of the user-
  1.1112 +    controllable parameters for colamd and symamd:
  1.1113 +
  1.1114 +        Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
  1.1115 +        entries are removed prior to ordering.  Columns with more than
  1.1116 +        max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
  1.1117 +        prior to ordering, and placed last in the output column ordering. 
  1.1118 +
  1.1119 +        Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
  1.1120 +        entries are removed prior to ordering, and placed last in the
  1.1121 +        output ordering.
  1.1122 +
  1.1123 +        knobs [0]       dense row control
  1.1124 +
  1.1125 +        knobs [1]       dense column control
  1.1126 +
  1.1127 +        knobs [2]       if nonzero, do aggresive absorption
  1.1128 +
  1.1129 +        knobs [3..19]   unused, but future versions might use this
  1.1130 +
  1.1131 +*/
  1.1132 +
  1.1133 +PUBLIC void COLAMD_set_defaults
  1.1134 +(
  1.1135 +    /* === Parameters ======================================================= */
  1.1136 +
  1.1137 +    double knobs [COLAMD_KNOBS]         /* knob array */
  1.1138 +)
  1.1139 +{
  1.1140 +    /* === Local variables ================================================== */
  1.1141 +
  1.1142 +    Int i ;
  1.1143 +
  1.1144 +    if (!knobs)
  1.1145 +    {
  1.1146 +        return ;                        /* no knobs to initialize */
  1.1147 +    }
  1.1148 +    for (i = 0 ; i < COLAMD_KNOBS ; i++)
  1.1149 +    {
  1.1150 +        knobs [i] = 0 ;
  1.1151 +    }
  1.1152 +    knobs [COLAMD_DENSE_ROW] = 10 ;
  1.1153 +    knobs [COLAMD_DENSE_COL] = 10 ;
  1.1154 +    knobs [COLAMD_AGGRESSIVE] = TRUE ;  /* default: do aggressive absorption*/
  1.1155 +}
  1.1156 +
  1.1157 +
  1.1158 +/* ========================================================================== */
  1.1159 +/* === symamd =============================================================== */
  1.1160 +/* ========================================================================== */
  1.1161 +
  1.1162 +PUBLIC Int SYMAMD_MAIN                  /* return TRUE if OK, FALSE otherwise */
  1.1163 +(
  1.1164 +    /* === Parameters ======================================================= */
  1.1165 +
  1.1166 +    Int n,                              /* number of rows and columns of A */
  1.1167 +    Int A [],                           /* row indices of A */
  1.1168 +    Int p [],                           /* column pointers of A */
  1.1169 +    Int perm [],                        /* output permutation, size n+1 */
  1.1170 +    double knobs [COLAMD_KNOBS],        /* parameters (uses defaults if NULL) */
  1.1171 +    Int stats [COLAMD_STATS],           /* output statistics and error codes */
  1.1172 +    void * (*allocate) (size_t, size_t),
  1.1173 +                                        /* pointer to calloc (ANSI C) or */
  1.1174 +                                        /* mxCalloc (for MATLAB mexFunction) */
  1.1175 +    void (*release) (void *)
  1.1176 +                                        /* pointer to free (ANSI C) or */
  1.1177 +                                        /* mxFree (for MATLAB mexFunction) */
  1.1178 +)
  1.1179 +{
  1.1180 +    /* === Local variables ================================================== */
  1.1181 +
  1.1182 +    Int *count ;                /* length of each column of M, and col pointer*/
  1.1183 +    Int *mark ;                 /* mark array for finding duplicate entries */
  1.1184 +    Int *M ;                    /* row indices of matrix M */
  1.1185 +    size_t Mlen ;               /* length of M */
  1.1186 +    Int n_row ;                 /* number of rows in M */
  1.1187 +    Int nnz ;                   /* number of entries in A */
  1.1188 +    Int i ;                     /* row index of A */
  1.1189 +    Int j ;                     /* column index of A */
  1.1190 +    Int k ;                     /* row index of M */ 
  1.1191 +    Int mnz ;                   /* number of nonzeros in M */
  1.1192 +    Int pp ;                    /* index into a column of A */
  1.1193 +    Int last_row ;              /* last row seen in the current column */
  1.1194 +    Int length ;                /* number of nonzeros in a column */
  1.1195 +
  1.1196 +    double cknobs [COLAMD_KNOBS] ;              /* knobs for colamd */
  1.1197 +    double default_knobs [COLAMD_KNOBS] ;       /* default knobs for colamd */
  1.1198 +
  1.1199 +#ifndef NDEBUG
  1.1200 +    colamd_get_debug ("symamd") ;
  1.1201 +#endif /* NDEBUG */
  1.1202 +
  1.1203 +    /* === Check the input arguments ======================================== */
  1.1204 +
  1.1205 +    if (!stats)
  1.1206 +    {
  1.1207 +        DEBUG0 (("symamd: stats not present\n")) ;
  1.1208 +        return (FALSE) ;
  1.1209 +    }
  1.1210 +    for (i = 0 ; i < COLAMD_STATS ; i++)
  1.1211 +    {
  1.1212 +        stats [i] = 0 ;
  1.1213 +    }
  1.1214 +    stats [COLAMD_STATUS] = COLAMD_OK ;
  1.1215 +    stats [COLAMD_INFO1] = -1 ;
  1.1216 +    stats [COLAMD_INFO2] = -1 ;
  1.1217 +
  1.1218 +    if (!A)
  1.1219 +    {
  1.1220 +        stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
  1.1221 +        DEBUG0 (("symamd: A not present\n")) ;
  1.1222 +        return (FALSE) ;
  1.1223 +    }
  1.1224 +
  1.1225 +    if (!p)             /* p is not present */
  1.1226 +    {
  1.1227 +        stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
  1.1228 +        DEBUG0 (("symamd: p not present\n")) ;
  1.1229 +        return (FALSE) ;
  1.1230 +    }
  1.1231 +
  1.1232 +    if (n < 0)          /* n must be >= 0 */
  1.1233 +    {
  1.1234 +        stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
  1.1235 +        stats [COLAMD_INFO1] = n ;
  1.1236 +        DEBUG0 (("symamd: n negative %d\n", n)) ;
  1.1237 +        return (FALSE) ;
  1.1238 +    }
  1.1239 +
  1.1240 +    nnz = p [n] ;
  1.1241 +    if (nnz < 0)        /* nnz must be >= 0 */
  1.1242 +    {
  1.1243 +        stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
  1.1244 +        stats [COLAMD_INFO1] = nnz ;
  1.1245 +        DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
  1.1246 +        return (FALSE) ;
  1.1247 +    }
  1.1248 +
  1.1249 +    if (p [0] != 0)
  1.1250 +    {
  1.1251 +        stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
  1.1252 +        stats [COLAMD_INFO1] = p [0] ;
  1.1253 +        DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
  1.1254 +        return (FALSE) ;
  1.1255 +    }
  1.1256 +
  1.1257 +    /* === If no knobs, set default knobs =================================== */
  1.1258 +
  1.1259 +    if (!knobs)
  1.1260 +    {
  1.1261 +        COLAMD_set_defaults (default_knobs) ;
  1.1262 +        knobs = default_knobs ;
  1.1263 +    }
  1.1264 +
  1.1265 +    /* === Allocate count and mark ========================================== */
  1.1266 +
  1.1267 +    count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
  1.1268 +    if (!count)
  1.1269 +    {
  1.1270 +        stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
  1.1271 +        DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
  1.1272 +        return (FALSE) ;
  1.1273 +    }
  1.1274 +
  1.1275 +    mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
  1.1276 +    if (!mark)
  1.1277 +    {
  1.1278 +        stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
  1.1279 +        (*release) ((void *) count) ;
  1.1280 +        DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
  1.1281 +        return (FALSE) ;
  1.1282 +    }
  1.1283 +
  1.1284 +    /* === Compute column counts of M, check if A is valid ================== */
  1.1285 +
  1.1286 +    stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
  1.1287 +
  1.1288 +    for (i = 0 ; i < n ; i++)
  1.1289 +    {
  1.1290 +        mark [i] = -1 ;
  1.1291 +    }
  1.1292 +
  1.1293 +    for (j = 0 ; j < n ; j++)
  1.1294 +    {
  1.1295 +        last_row = -1 ;
  1.1296 +
  1.1297 +        length = p [j+1] - p [j] ;
  1.1298 +        if (length < 0)
  1.1299 +        {
  1.1300 +            /* column pointers must be non-decreasing */
  1.1301 +            stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
  1.1302 +            stats [COLAMD_INFO1] = j ;
  1.1303 +            stats [COLAMD_INFO2] = length ;
  1.1304 +            (*release) ((void *) count) ;
  1.1305 +            (*release) ((void *) mark) ;
  1.1306 +            DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
  1.1307 +            return (FALSE) ;
  1.1308 +        }
  1.1309 +
  1.1310 +        for (pp = p [j] ; pp < p [j+1] ; pp++)
  1.1311 +        {
  1.1312 +            i = A [pp] ;
  1.1313 +            if (i < 0 || i >= n)
  1.1314 +            {
  1.1315 +                /* row index i, in column j, is out of bounds */
  1.1316 +                stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
  1.1317 +                stats [COLAMD_INFO1] = j ;
  1.1318 +                stats [COLAMD_INFO2] = i ;
  1.1319 +                stats [COLAMD_INFO3] = n ;
  1.1320 +                (*release) ((void *) count) ;
  1.1321 +                (*release) ((void *) mark) ;
  1.1322 +                DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
  1.1323 +                return (FALSE) ;
  1.1324 +            }
  1.1325 +
  1.1326 +            if (i <= last_row || mark [i] == j)
  1.1327 +            {
  1.1328 +                /* row index is unsorted or repeated (or both), thus col */
  1.1329 +                /* is jumbled.  This is a notice, not an error condition. */
  1.1330 +                stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
  1.1331 +                stats [COLAMD_INFO1] = j ;
  1.1332 +                stats [COLAMD_INFO2] = i ;
  1.1333 +                (stats [COLAMD_INFO3]) ++ ;
  1.1334 +                DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
  1.1335 +            }
  1.1336 +
  1.1337 +            if (i > j && mark [i] != j)
  1.1338 +            {
  1.1339 +                /* row k of M will contain column indices i and j */
  1.1340 +                count [i]++ ;
  1.1341 +                count [j]++ ;
  1.1342 +            }
  1.1343 +
  1.1344 +            /* mark the row as having been seen in this column */
  1.1345 +            mark [i] = j ;
  1.1346 +
  1.1347 +            last_row = i ;
  1.1348 +        }
  1.1349 +    }
  1.1350 +
  1.1351 +    /* v2.4: removed free(mark) */
  1.1352 +
  1.1353 +    /* === Compute column pointers of M ===================================== */
  1.1354 +
  1.1355 +    /* use output permutation, perm, for column pointers of M */
  1.1356 +    perm [0] = 0 ;
  1.1357 +    for (j = 1 ; j <= n ; j++)
  1.1358 +    {
  1.1359 +        perm [j] = perm [j-1] + count [j-1] ;
  1.1360 +    }
  1.1361 +    for (j = 0 ; j < n ; j++)
  1.1362 +    {
  1.1363 +        count [j] = perm [j] ;
  1.1364 +    }
  1.1365 +
  1.1366 +    /* === Construct M ====================================================== */
  1.1367 +
  1.1368 +    mnz = perm [n] ;
  1.1369 +    n_row = mnz / 2 ;
  1.1370 +    Mlen = COLAMD_recommended (mnz, n_row, n) ;
  1.1371 +    M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
  1.1372 +    DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
  1.1373 +        n_row, n, mnz, (double) Mlen)) ;
  1.1374 +
  1.1375 +    if (!M)
  1.1376 +    {
  1.1377 +        stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
  1.1378 +        (*release) ((void *) count) ;
  1.1379 +        (*release) ((void *) mark) ;
  1.1380 +        DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
  1.1381 +        return (FALSE) ;
  1.1382 +    }
  1.1383 +
  1.1384 +    k = 0 ;
  1.1385 +
  1.1386 +    if (stats [COLAMD_STATUS] == COLAMD_OK)
  1.1387 +    {
  1.1388 +        /* Matrix is OK */
  1.1389 +        for (j = 0 ; j < n ; j++)
  1.1390 +        {
  1.1391 +            ASSERT (p [j+1] - p [j] >= 0) ;
  1.1392 +            for (pp = p [j] ; pp < p [j+1] ; pp++)
  1.1393 +            {
  1.1394 +                i = A [pp] ;
  1.1395 +                ASSERT (i >= 0 && i < n) ;
  1.1396 +                if (i > j)
  1.1397 +                {
  1.1398 +                    /* row k of M contains column indices i and j */
  1.1399 +                    M [count [i]++] = k ;
  1.1400 +                    M [count [j]++] = k ;
  1.1401 +                    k++ ;
  1.1402 +                }
  1.1403 +            }
  1.1404 +        }
  1.1405 +    }
  1.1406 +    else
  1.1407 +    {
  1.1408 +        /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
  1.1409 +        DEBUG0 (("symamd: Duplicates in A.\n")) ;
  1.1410 +        for (i = 0 ; i < n ; i++)
  1.1411 +        {
  1.1412 +            mark [i] = -1 ;
  1.1413 +        }
  1.1414 +        for (j = 0 ; j < n ; j++)
  1.1415 +        {
  1.1416 +            ASSERT (p [j+1] - p [j] >= 0) ;
  1.1417 +            for (pp = p [j] ; pp < p [j+1] ; pp++)
  1.1418 +            {
  1.1419 +                i = A [pp] ;
  1.1420 +                ASSERT (i >= 0 && i < n) ;
  1.1421 +                if (i > j && mark [i] != j)
  1.1422 +                {
  1.1423 +                    /* row k of M contains column indices i and j */
  1.1424 +                    M [count [i]++] = k ;
  1.1425 +                    M [count [j]++] = k ;
  1.1426 +                    k++ ;
  1.1427 +                    mark [i] = j ;
  1.1428 +                }
  1.1429 +            }
  1.1430 +        }
  1.1431 +        /* v2.4: free(mark) moved below */
  1.1432 +    }
  1.1433 +
  1.1434 +    /* count and mark no longer needed */
  1.1435 +    (*release) ((void *) count) ;
  1.1436 +    (*release) ((void *) mark) ;        /* v2.4: free (mark) moved here */
  1.1437 +    ASSERT (k == n_row) ;
  1.1438 +
  1.1439 +    /* === Adjust the knobs for M =========================================== */
  1.1440 +
  1.1441 +    for (i = 0 ; i < COLAMD_KNOBS ; i++)
  1.1442 +    {
  1.1443 +        cknobs [i] = knobs [i] ;
  1.1444 +    }
  1.1445 +
  1.1446 +    /* there are no dense rows in M */
  1.1447 +    cknobs [COLAMD_DENSE_ROW] = -1 ;
  1.1448 +    cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
  1.1449 +
  1.1450 +    /* === Order the columns of M =========================================== */
  1.1451 +
  1.1452 +    /* v2.4: colamd cannot fail here, so the error check is removed */
  1.1453 +    (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
  1.1454 +
  1.1455 +    /* Note that the output permutation is now in perm */
  1.1456 +
  1.1457 +    /* === get the statistics for symamd from colamd ======================== */
  1.1458 +
  1.1459 +    /* a dense column in colamd means a dense row and col in symamd */
  1.1460 +    stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
  1.1461 +
  1.1462 +    /* === Free M =========================================================== */
  1.1463 +
  1.1464 +    (*release) ((void *) M) ;
  1.1465 +    DEBUG0 (("symamd: done.\n")) ;
  1.1466 +    return (TRUE) ;
  1.1467 +
  1.1468 +}
  1.1469 +
  1.1470 +/* ========================================================================== */
  1.1471 +/* === colamd =============================================================== */
  1.1472 +/* ========================================================================== */
  1.1473 +
  1.1474 +/*
  1.1475 +    The colamd routine computes a column ordering Q of a sparse matrix
  1.1476 +    A such that the LU factorization P(AQ) = LU remains sparse, where P is
  1.1477 +    selected via partial pivoting.   The routine can also be viewed as
  1.1478 +    providing a permutation Q such that the Cholesky factorization
  1.1479 +    (AQ)'(AQ) = LL' remains sparse.
  1.1480 +*/
  1.1481 +
  1.1482 +PUBLIC Int COLAMD_MAIN          /* returns TRUE if successful, FALSE otherwise*/
  1.1483 +(
  1.1484 +    /* === Parameters ======================================================= */
  1.1485 +
  1.1486 +    Int n_row,                  /* number of rows in A */
  1.1487 +    Int n_col,                  /* number of columns in A */
  1.1488 +    Int Alen,                   /* length of A */
  1.1489 +    Int A [],                   /* row indices of A */
  1.1490 +    Int p [],                   /* pointers to columns in A */
  1.1491 +    double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
  1.1492 +    Int stats [COLAMD_STATS]    /* output statistics and error codes */
  1.1493 +)
  1.1494 +{
  1.1495 +    /* === Local variables ================================================== */
  1.1496 +
  1.1497 +    Int i ;                     /* loop index */
  1.1498 +    Int nnz ;                   /* nonzeros in A */
  1.1499 +    size_t Row_size ;           /* size of Row [], in integers */
  1.1500 +    size_t Col_size ;           /* size of Col [], in integers */
  1.1501 +    size_t need ;               /* minimum required length of A */
  1.1502 +    Colamd_Row *Row ;           /* pointer into A of Row [0..n_row] array */
  1.1503 +    Colamd_Col *Col ;           /* pointer into A of Col [0..n_col] array */
  1.1504 +    Int n_col2 ;                /* number of non-dense, non-empty columns */
  1.1505 +    Int n_row2 ;                /* number of non-dense, non-empty rows */
  1.1506 +    Int ngarbage ;              /* number of garbage collections performed */
  1.1507 +    Int max_deg ;               /* maximum row degree */
  1.1508 +    double default_knobs [COLAMD_KNOBS] ;       /* default knobs array */
  1.1509 +    Int aggressive ;            /* do aggressive absorption */
  1.1510 +    int ok ;
  1.1511 +
  1.1512 +#ifndef NDEBUG
  1.1513 +    colamd_get_debug ("colamd") ;
  1.1514 +#endif /* NDEBUG */
  1.1515 +
  1.1516 +    /* === Check the input arguments ======================================== */
  1.1517 +
  1.1518 +    if (!stats)
  1.1519 +    {
  1.1520 +        DEBUG0 (("colamd: stats not present\n")) ;
  1.1521 +        return (FALSE) ;
  1.1522 +    }
  1.1523 +    for (i = 0 ; i < COLAMD_STATS ; i++)
  1.1524 +    {
  1.1525 +        stats [i] = 0 ;
  1.1526 +    }
  1.1527 +    stats [COLAMD_STATUS] = COLAMD_OK ;
  1.1528 +    stats [COLAMD_INFO1] = -1 ;
  1.1529 +    stats [COLAMD_INFO2] = -1 ;
  1.1530 +
  1.1531 +    if (!A)             /* A is not present */
  1.1532 +    {
  1.1533 +        stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
  1.1534 +        DEBUG0 (("colamd: A not present\n")) ;
  1.1535 +        return (FALSE) ;
  1.1536 +    }
  1.1537 +
  1.1538 +    if (!p)             /* p is not present */
  1.1539 +    {
  1.1540 +        stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
  1.1541 +        DEBUG0 (("colamd: p not present\n")) ;
  1.1542 +        return (FALSE) ;
  1.1543 +    }
  1.1544 +
  1.1545 +    if (n_row < 0)      /* n_row must be >= 0 */
  1.1546 +    {
  1.1547 +        stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
  1.1548 +        stats [COLAMD_INFO1] = n_row ;
  1.1549 +        DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
  1.1550 +        return (FALSE) ;
  1.1551 +    }
  1.1552 +
  1.1553 +    if (n_col < 0)      /* n_col must be >= 0 */
  1.1554 +    {
  1.1555 +        stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
  1.1556 +        stats [COLAMD_INFO1] = n_col ;
  1.1557 +        DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
  1.1558 +        return (FALSE) ;
  1.1559 +    }
  1.1560 +
  1.1561 +    nnz = p [n_col] ;
  1.1562 +    if (nnz < 0)        /* nnz must be >= 0 */
  1.1563 +    {
  1.1564 +        stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
  1.1565 +        stats [COLAMD_INFO1] = nnz ;
  1.1566 +        DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
  1.1567 +        return (FALSE) ;
  1.1568 +    }
  1.1569 +
  1.1570 +    if (p [0] != 0)
  1.1571 +    {
  1.1572 +        stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
  1.1573 +        stats [COLAMD_INFO1] = p [0] ;
  1.1574 +        DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
  1.1575 +        return (FALSE) ;
  1.1576 +    }
  1.1577 +
  1.1578 +    /* === If no knobs, set default knobs =================================== */
  1.1579 +
  1.1580 +    if (!knobs)
  1.1581 +    {
  1.1582 +        COLAMD_set_defaults (default_knobs) ;
  1.1583 +        knobs = default_knobs ;
  1.1584 +    }
  1.1585 +
  1.1586 +    aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
  1.1587 +
  1.1588 +    /* === Allocate the Row and Col arrays from array A ===================== */
  1.1589 +
  1.1590 +    ok = TRUE ;
  1.1591 +    Col_size = COLAMD_C (n_col, &ok) ;      /* size of Col array of structs */
  1.1592 +    Row_size = COLAMD_R (n_row, &ok) ;      /* size of Row array of structs */
  1.1593 +
  1.1594 +    /* need = 2*nnz + n_col + Col_size + Row_size ; */
  1.1595 +    need = t_mult (nnz, 2, &ok) ;
  1.1596 +    need = t_add (need, n_col, &ok) ;
  1.1597 +    need = t_add (need, Col_size, &ok) ;
  1.1598 +    need = t_add (need, Row_size, &ok) ;
  1.1599 +
  1.1600 +    if (!ok || need > (size_t) Alen || need > Int_MAX)
  1.1601 +    {
  1.1602 +        /* not enough space in array A to perform the ordering */
  1.1603 +        stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
  1.1604 +        stats [COLAMD_INFO1] = need ;
  1.1605 +        stats [COLAMD_INFO2] = Alen ;
  1.1606 +        DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
  1.1607 +        return (FALSE) ;
  1.1608 +    }
  1.1609 +
  1.1610 +    Alen -= Col_size + Row_size ;
  1.1611 +    Col = (Colamd_Col *) &A [Alen] ;
  1.1612 +    Row = (Colamd_Row *) &A [Alen + Col_size] ;
  1.1613 +
  1.1614 +    /* === Construct the row and column data structures ===================== */
  1.1615 +
  1.1616 +    if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
  1.1617 +    {
  1.1618 +        /* input matrix is invalid */
  1.1619 +        DEBUG0 (("colamd: Matrix invalid\n")) ;
  1.1620 +        return (FALSE) ;
  1.1621 +    }
  1.1622 +
  1.1623 +    /* === Initialize scores, kill dense rows/columns ======================= */
  1.1624 +
  1.1625 +    init_scoring (n_row, n_col, Row, Col, A, p, knobs,
  1.1626 +        &n_row2, &n_col2, &max_deg) ;
  1.1627 +
  1.1628 +    /* === Order the supercolumns =========================================== */
  1.1629 +
  1.1630 +    ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
  1.1631 +        n_col2, max_deg, 2*nnz, aggressive) ;
  1.1632 +
  1.1633 +    /* === Order the non-principal columns ================================== */
  1.1634 +
  1.1635 +    order_children (n_col, Col, p) ;
  1.1636 +
  1.1637 +    /* === Return statistics in stats ======================================= */
  1.1638 +
  1.1639 +    stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
  1.1640 +    stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
  1.1641 +    stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
  1.1642 +    DEBUG0 (("colamd: done.\n")) ; 
  1.1643 +    return (TRUE) ;
  1.1644 +}
  1.1645 +
  1.1646 +
  1.1647 +/* ========================================================================== */
  1.1648 +/* === colamd_report ======================================================== */
  1.1649 +/* ========================================================================== */
  1.1650 +
  1.1651 +PUBLIC void COLAMD_report
  1.1652 +(
  1.1653 +    Int stats [COLAMD_STATS]
  1.1654 +)
  1.1655 +{
  1.1656 +    print_report ("colamd", stats) ;
  1.1657 +}
  1.1658 +
  1.1659 +
  1.1660 +/* ========================================================================== */
  1.1661 +/* === symamd_report ======================================================== */
  1.1662 +/* ========================================================================== */
  1.1663 +
  1.1664 +PUBLIC void SYMAMD_report
  1.1665 +(
  1.1666 +    Int stats [COLAMD_STATS]
  1.1667 +)
  1.1668 +{
  1.1669 +    print_report ("symamd", stats) ;
  1.1670 +}
  1.1671 +
  1.1672 +
  1.1673 +
  1.1674 +/* ========================================================================== */
  1.1675 +/* === NON-USER-CALLABLE ROUTINES: ========================================== */
  1.1676 +/* ========================================================================== */
  1.1677 +
  1.1678 +/* There are no user-callable routines beyond this point in the file */
  1.1679 +
  1.1680 +
  1.1681 +/* ========================================================================== */
  1.1682 +/* === init_rows_cols ======================================================= */
  1.1683 +/* ========================================================================== */
  1.1684 +
  1.1685 +/*
  1.1686 +    Takes the column form of the matrix in A and creates the row form of the
  1.1687 +    matrix.  Also, row and column attributes are stored in the Col and Row
  1.1688 +    structs.  If the columns are un-sorted or contain duplicate row indices,
  1.1689 +    this routine will also sort and remove duplicate row indices from the
  1.1690 +    column form of the matrix.  Returns FALSE if the matrix is invalid,
  1.1691 +    TRUE otherwise.  Not user-callable.
  1.1692 +*/
  1.1693 +
  1.1694 +PRIVATE Int init_rows_cols      /* returns TRUE if OK, or FALSE otherwise */
  1.1695 +(
  1.1696 +    /* === Parameters ======================================================= */
  1.1697 +
  1.1698 +    Int n_row,                  /* number of rows of A */
  1.1699 +    Int n_col,                  /* number of columns of A */
  1.1700 +    Colamd_Row Row [],          /* of size n_row+1 */
  1.1701 +    Colamd_Col Col [],          /* of size n_col+1 */
  1.1702 +    Int A [],                   /* row indices of A, of size Alen */
  1.1703 +    Int p [],                   /* pointers to columns in A, of size n_col+1 */
  1.1704 +    Int stats [COLAMD_STATS]    /* colamd statistics */ 
  1.1705 +)
  1.1706 +{
  1.1707 +    /* === Local variables ================================================== */
  1.1708 +
  1.1709 +    Int col ;                   /* a column index */
  1.1710 +    Int row ;                   /* a row index */
  1.1711 +    Int *cp ;                   /* a column pointer */
  1.1712 +    Int *cp_end ;               /* a pointer to the end of a column */
  1.1713 +    Int *rp ;                   /* a row pointer */
  1.1714 +    Int *rp_end ;               /* a pointer to the end of a row */
  1.1715 +    Int last_row ;              /* previous row */
  1.1716 +
  1.1717 +    /* === Initialize columns, and check column pointers ==================== */
  1.1718 +
  1.1719 +    for (col = 0 ; col < n_col ; col++)
  1.1720 +    {
  1.1721 +        Col [col].start = p [col] ;
  1.1722 +        Col [col].length = p [col+1] - p [col] ;
  1.1723 +
  1.1724 +        if (Col [col].length < 0)
  1.1725 +        {
  1.1726 +            /* column pointers must be non-decreasing */
  1.1727 +            stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
  1.1728 +            stats [COLAMD_INFO1] = col ;
  1.1729 +            stats [COLAMD_INFO2] = Col [col].length ;
  1.1730 +            DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
  1.1731 +            return (FALSE) ;
  1.1732 +        }
  1.1733 +
  1.1734 +        Col [col].shared1.thickness = 1 ;
  1.1735 +        Col [col].shared2.score = 0 ;
  1.1736 +        Col [col].shared3.prev = EMPTY ;
  1.1737 +        Col [col].shared4.degree_next = EMPTY ;
  1.1738 +    }
  1.1739 +
  1.1740 +    /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
  1.1741 +
  1.1742 +    /* === Scan columns, compute row degrees, and check row indices ========= */
  1.1743 +
  1.1744 +    stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
  1.1745 +
  1.1746 +    for (row = 0 ; row < n_row ; row++)
  1.1747 +    {
  1.1748 +        Row [row].length = 0 ;
  1.1749 +        Row [row].shared2.mark = -1 ;
  1.1750 +    }
  1.1751 +
  1.1752 +    for (col = 0 ; col < n_col ; col++)
  1.1753 +    {
  1.1754 +        last_row = -1 ;
  1.1755 +
  1.1756 +        cp = &A [p [col]] ;
  1.1757 +        cp_end = &A [p [col+1]] ;
  1.1758 +
  1.1759 +        while (cp < cp_end)
  1.1760 +        {
  1.1761 +            row = *cp++ ;
  1.1762 +
  1.1763 +            /* make sure row indices within range */
  1.1764 +            if (row < 0 || row >= n_row)
  1.1765 +            {
  1.1766 +                stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
  1.1767 +                stats [COLAMD_INFO1] = col ;
  1.1768 +                stats [COLAMD_INFO2] = row ;
  1.1769 +                stats [COLAMD_INFO3] = n_row ;
  1.1770 +                DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
  1.1771 +                return (FALSE) ;
  1.1772 +            }
  1.1773 +
  1.1774 +            if (row <= last_row || Row [row].shared2.mark == col)
  1.1775 +            {
  1.1776 +                /* row index are unsorted or repeated (or both), thus col */
  1.1777 +                /* is jumbled.  This is a notice, not an error condition. */
  1.1778 +                stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
  1.1779 +                stats [COLAMD_INFO1] = col ;
  1.1780 +                stats [COLAMD_INFO2] = row ;
  1.1781 +                (stats [COLAMD_INFO3]) ++ ;
  1.1782 +                DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
  1.1783 +            }
  1.1784 +
  1.1785 +            if (Row [row].shared2.mark != col)
  1.1786 +            {
  1.1787 +                Row [row].length++ ;
  1.1788 +            }
  1.1789 +            else
  1.1790 +            {
  1.1791 +                /* this is a repeated entry in the column, */
  1.1792 +                /* it will be removed */
  1.1793 +                Col [col].length-- ;
  1.1794 +            }
  1.1795 +
  1.1796 +            /* mark the row as having been seen in this column */
  1.1797 +            Row [row].shared2.mark = col ;
  1.1798 +
  1.1799 +            last_row = row ;
  1.1800 +        }
  1.1801 +    }
  1.1802 +
  1.1803 +    /* === Compute row pointers ============================================= */
  1.1804 +
  1.1805 +    /* row form of the matrix starts directly after the column */
  1.1806 +    /* form of matrix in A */
  1.1807 +    Row [0].start = p [n_col] ;
  1.1808 +    Row [0].shared1.p = Row [0].start ;
  1.1809 +    Row [0].shared2.mark = -1 ;
  1.1810 +    for (row = 1 ; row < n_row ; row++)
  1.1811 +    {
  1.1812 +        Row [row].start = Row [row-1].start + Row [row-1].length ;
  1.1813 +        Row [row].shared1.p = Row [row].start ;
  1.1814 +        Row [row].shared2.mark = -1 ;
  1.1815 +    }
  1.1816 +
  1.1817 +    /* === Create row form ================================================== */
  1.1818 +
  1.1819 +    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
  1.1820 +    {
  1.1821 +        /* if cols jumbled, watch for repeated row indices */
  1.1822 +        for (col = 0 ; col < n_col ; col++)
  1.1823 +        {
  1.1824 +            cp = &A [p [col]] ;
  1.1825 +            cp_end = &A [p [col+1]] ;
  1.1826 +            while (cp < cp_end)
  1.1827 +            {
  1.1828 +                row = *cp++ ;
  1.1829 +                if (Row [row].shared2.mark != col)
  1.1830 +                {
  1.1831 +                    A [(Row [row].shared1.p)++] = col ;
  1.1832 +                    Row [row].shared2.mark = col ;
  1.1833 +                }
  1.1834 +            }
  1.1835 +        }
  1.1836 +    }
  1.1837 +    else
  1.1838 +    {
  1.1839 +        /* if cols not jumbled, we don't need the mark (this is faster) */
  1.1840 +        for (col = 0 ; col < n_col ; col++)
  1.1841 +        {
  1.1842 +            cp = &A [p [col]] ;
  1.1843 +            cp_end = &A [p [col+1]] ;
  1.1844 +            while (cp < cp_end)
  1.1845 +            {
  1.1846 +                A [(Row [*cp++].shared1.p)++] = col ;
  1.1847 +            }
  1.1848 +        }
  1.1849 +    }
  1.1850 +
  1.1851 +    /* === Clear the row marks and set row degrees ========================== */
  1.1852 +
  1.1853 +    for (row = 0 ; row < n_row ; row++)
  1.1854 +    {
  1.1855 +        Row [row].shared2.mark = 0 ;
  1.1856 +        Row [row].shared1.degree = Row [row].length ;
  1.1857 +    }
  1.1858 +
  1.1859 +    /* === See if we need to re-create columns ============================== */
  1.1860 +
  1.1861 +    if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
  1.1862 +    {
  1.1863 +        DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
  1.1864 +
  1.1865 +#ifndef NDEBUG
  1.1866 +        /* make sure column lengths are correct */
  1.1867 +        for (col = 0 ; col < n_col ; col++)
  1.1868 +        {
  1.1869 +            p [col] = Col [col].length ;
  1.1870 +        }
  1.1871 +        for (row = 0 ; row < n_row ; row++)
  1.1872 +        {
  1.1873 +            rp = &A [Row [row].start] ;
  1.1874 +            rp_end = rp + Row [row].length ;
  1.1875 +            while (rp < rp_end)
  1.1876 +            {
  1.1877 +                p [*rp++]-- ;
  1.1878 +            }
  1.1879 +        }
  1.1880 +        for (col = 0 ; col < n_col ; col++)
  1.1881 +        {
  1.1882 +            ASSERT (p [col] == 0) ;
  1.1883 +        }
  1.1884 +        /* now p is all zero (different than when debugging is turned off) */
  1.1885 +#endif /* NDEBUG */
  1.1886 +
  1.1887 +        /* === Compute col pointers ========================================= */
  1.1888 +
  1.1889 +        /* col form of the matrix starts at A [0]. */
  1.1890 +        /* Note, we may have a gap between the col form and the row */
  1.1891 +        /* form if there were duplicate entries, if so, it will be */
  1.1892 +        /* removed upon the first garbage collection */
  1.1893 +        Col [0].start = 0 ;
  1.1894 +        p [0] = Col [0].start ;
  1.1895 +        for (col = 1 ; col < n_col ; col++)
  1.1896 +        {
  1.1897 +            /* note that the lengths here are for pruned columns, i.e. */
  1.1898 +            /* no duplicate row indices will exist for these columns */
  1.1899 +            Col [col].start = Col [col-1].start + Col [col-1].length ;
  1.1900 +            p [col] = Col [col].start ;
  1.1901 +        }
  1.1902 +
  1.1903 +        /* === Re-create col form =========================================== */
  1.1904 +
  1.1905 +        for (row = 0 ; row < n_row ; row++)
  1.1906 +        {
  1.1907 +            rp = &A [Row [row].start] ;
  1.1908 +            rp_end = rp + Row [row].length ;
  1.1909 +            while (rp < rp_end)
  1.1910 +            {
  1.1911 +                A [(p [*rp++])++] = row ;
  1.1912 +            }
  1.1913 +        }
  1.1914 +    }
  1.1915 +
  1.1916 +    /* === Done.  Matrix is not (or no longer) jumbled ====================== */
  1.1917 +
  1.1918 +    return (TRUE) ;
  1.1919 +}
  1.1920 +
  1.1921 +
  1.1922 +/* ========================================================================== */
  1.1923 +/* === init_scoring ========================================================= */
  1.1924 +/* ========================================================================== */
  1.1925 +
  1.1926 +/*
  1.1927 +    Kills dense or empty columns and rows, calculates an initial score for
  1.1928 +    each column, and places all columns in the degree lists.  Not user-callable.
  1.1929 +*/
  1.1930 +
  1.1931 +PRIVATE void init_scoring
  1.1932 +(
  1.1933 +    /* === Parameters ======================================================= */
  1.1934 +
  1.1935 +    Int n_row,                  /* number of rows of A */
  1.1936 +    Int n_col,                  /* number of columns of A */
  1.1937 +    Colamd_Row Row [],          /* of size n_row+1 */
  1.1938 +    Colamd_Col Col [],          /* of size n_col+1 */
  1.1939 +    Int A [],                   /* column form and row form of A */
  1.1940 +    Int head [],                /* of size n_col+1 */
  1.1941 +    double knobs [COLAMD_KNOBS],/* parameters */
  1.1942 +    Int *p_n_row2,              /* number of non-dense, non-empty rows */
  1.1943 +    Int *p_n_col2,              /* number of non-dense, non-empty columns */
  1.1944 +    Int *p_max_deg              /* maximum row degree */
  1.1945 +)
  1.1946 +{
  1.1947 +    /* === Local variables ================================================== */
  1.1948 +
  1.1949 +    Int c ;                     /* a column index */
  1.1950 +    Int r, row ;                /* a row index */
  1.1951 +    Int *cp ;                   /* a column pointer */
  1.1952 +    Int deg ;                   /* degree of a row or column */
  1.1953 +    Int *cp_end ;               /* a pointer to the end of a column */
  1.1954 +    Int *new_cp ;               /* new column pointer */
  1.1955 +    Int col_length ;            /* length of pruned column */
  1.1956 +    Int score ;                 /* current column score */
  1.1957 +    Int n_col2 ;                /* number of non-dense, non-empty columns */
  1.1958 +    Int n_row2 ;                /* number of non-dense, non-empty rows */
  1.1959 +    Int dense_row_count ;       /* remove rows with more entries than this */
  1.1960 +    Int dense_col_count ;       /* remove cols with more entries than this */
  1.1961 +    Int min_score ;             /* smallest column score */
  1.1962 +    Int max_deg ;               /* maximum row degree */
  1.1963 +    Int next_col ;              /* Used to add to degree list.*/
  1.1964 +
  1.1965 +#ifndef NDEBUG
  1.1966 +    Int debug_count ;           /* debug only. */
  1.1967 +#endif /* NDEBUG */
  1.1968 +
  1.1969 +    /* === Extract knobs ==================================================== */
  1.1970 +
  1.1971 +    /* Note: if knobs contains a NaN, this is undefined: */
  1.1972 +    if (knobs [COLAMD_DENSE_ROW] < 0)
  1.1973 +    {
  1.1974 +        /* only remove completely dense rows */
  1.1975 +        dense_row_count = n_col-1 ;
  1.1976 +    }
  1.1977 +    else
  1.1978 +    {
  1.1979 +        dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
  1.1980 +    }
  1.1981 +    if (knobs [COLAMD_DENSE_COL] < 0)
  1.1982 +    {
  1.1983 +        /* only remove completely dense columns */
  1.1984 +        dense_col_count = n_row-1 ;
  1.1985 +    }
  1.1986 +    else
  1.1987 +    {
  1.1988 +        dense_col_count =
  1.1989 +            DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
  1.1990 +    }
  1.1991 +
  1.1992 +    DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
  1.1993 +    max_deg = 0 ;
  1.1994 +    n_col2 = n_col ;
  1.1995 +    n_row2 = n_row ;
  1.1996 +
  1.1997 +    /* === Kill empty columns =============================================== */
  1.1998 +
  1.1999 +    /* Put the empty columns at the end in their natural order, so that LU */
  1.2000 +    /* factorization can proceed as far as possible. */
  1.2001 +    for (c = n_col-1 ; c >= 0 ; c--)
  1.2002 +    {
  1.2003 +        deg = Col [c].length ;
  1.2004 +        if (deg == 0)
  1.2005 +        {
  1.2006 +            /* this is a empty column, kill and order it last */
  1.2007 +            Col [c].shared2.order = --n_col2 ;
  1.2008 +            KILL_PRINCIPAL_COL (c) ;
  1.2009 +        }
  1.2010 +    }
  1.2011 +    DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
  1.2012 +
  1.2013 +    /* === Kill dense columns =============================================== */
  1.2014 +
  1.2015 +    /* Put the dense columns at the end, in their natural order */
  1.2016 +    for (c = n_col-1 ; c >= 0 ; c--)
  1.2017 +    {
  1.2018 +        /* skip any dead columns */
  1.2019 +        if (COL_IS_DEAD (c))
  1.2020 +        {
  1.2021 +            continue ;
  1.2022 +        }
  1.2023 +        deg = Col [c].length ;
  1.2024 +        if (deg > dense_col_count)
  1.2025 +        {
  1.2026 +            /* this is a dense column, kill and order it last */
  1.2027 +            Col [c].shared2.order = --n_col2 ;
  1.2028 +            /* decrement the row degrees */
  1.2029 +            cp = &A [Col [c].start] ;
  1.2030 +            cp_end = cp + Col [c].length ;
  1.2031 +            while (cp < cp_end)
  1.2032 +            {
  1.2033 +                Row [*cp++].shared1.degree-- ;
  1.2034 +            }
  1.2035 +            KILL_PRINCIPAL_COL (c) ;
  1.2036 +        }
  1.2037 +    }
  1.2038 +    DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
  1.2039 +
  1.2040 +    /* === Kill dense and empty rows ======================================== */
  1.2041 +
  1.2042 +    for (r = 0 ; r < n_row ; r++)
  1.2043 +    {
  1.2044 +        deg = Row [r].shared1.degree ;
  1.2045 +        ASSERT (deg >= 0 && deg <= n_col) ;
  1.2046 +        if (deg > dense_row_count || deg == 0)
  1.2047 +        {
  1.2048 +            /* kill a dense or empty row */
  1.2049 +            KILL_ROW (r) ;
  1.2050 +            --n_row2 ;
  1.2051 +        }
  1.2052 +        else
  1.2053 +        {
  1.2054 +            /* keep track of max degree of remaining rows */
  1.2055 +            max_deg = MAX (max_deg, deg) ;
  1.2056 +        }
  1.2057 +    }
  1.2058 +    DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
  1.2059 +
  1.2060 +    /* === Compute initial column scores ==================================== */
  1.2061 +
  1.2062 +    /* At this point the row degrees are accurate.  They reflect the number */
  1.2063 +    /* of "live" (non-dense) columns in each row.  No empty rows exist. */
  1.2064 +    /* Some "live" columns may contain only dead rows, however.  These are */
  1.2065 +    /* pruned in the code below. */
  1.2066 +
  1.2067 +    /* now find the initial matlab score for each column */
  1.2068 +    for (c = n_col-1 ; c >= 0 ; c--)
  1.2069 +    {
  1.2070 +        /* skip dead column */
  1.2071 +        if (COL_IS_DEAD (c))
  1.2072 +        {
  1.2073 +            continue ;
  1.2074 +        }
  1.2075 +        score = 0 ;
  1.2076 +        cp = &A [Col [c].start] ;
  1.2077 +        new_cp = cp ;
  1.2078 +        cp_end = cp + Col [c].length ;
  1.2079 +        while (cp < cp_end)
  1.2080 +        {
  1.2081 +            /* get a row */
  1.2082 +            row = *cp++ ;
  1.2083 +            /* skip if dead */
  1.2084 +            if (ROW_IS_DEAD (row))
  1.2085 +            {
  1.2086 +                continue ;
  1.2087 +            }
  1.2088 +            /* compact the column */
  1.2089 +            *new_cp++ = row ;
  1.2090 +            /* add row's external degree */
  1.2091 +            score += Row [row].shared1.degree - 1 ;
  1.2092 +            /* guard against integer overflow */
  1.2093 +            score = MIN (score, n_col) ;
  1.2094 +        }
  1.2095 +        /* determine pruned column length */
  1.2096 +        col_length = (Int) (new_cp - &A [Col [c].start]) ;
  1.2097 +        if (col_length == 0)
  1.2098 +        {
  1.2099 +            /* a newly-made null column (all rows in this col are "dense" */
  1.2100 +            /* and have already been killed) */
  1.2101 +            DEBUG2 (("Newly null killed: %d\n", c)) ;
  1.2102 +            Col [c].shared2.order = --n_col2 ;
  1.2103 +            KILL_PRINCIPAL_COL (c) ;
  1.2104 +        }
  1.2105 +        else
  1.2106 +        {
  1.2107 +            /* set column length and set score */
  1.2108 +            ASSERT (score >= 0) ;
  1.2109 +            ASSERT (score <= n_col) ;
  1.2110 +            Col [c].length = col_length ;
  1.2111 +            Col [c].shared2.score = score ;
  1.2112 +        }
  1.2113 +    }
  1.2114 +    DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
  1.2115 +        n_col-n_col2)) ;
  1.2116 +
  1.2117 +    /* At this point, all empty rows and columns are dead.  All live columns */
  1.2118 +    /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
  1.2119 +    /* yet).  Rows may contain dead columns, but all live rows contain at */
  1.2120 +    /* least one live column. */
  1.2121 +
  1.2122 +#ifndef NDEBUG
  1.2123 +    debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
  1.2124 +#endif /* NDEBUG */
  1.2125 +
  1.2126 +    /* === Initialize degree lists ========================================== */
  1.2127 +
  1.2128 +#ifndef NDEBUG
  1.2129 +    debug_count = 0 ;
  1.2130 +#endif /* NDEBUG */
  1.2131 +
  1.2132 +    /* clear the hash buckets */
  1.2133 +    for (c = 0 ; c <= n_col ; c++)
  1.2134 +    {
  1.2135 +        head [c] = EMPTY ;
  1.2136 +    }
  1.2137 +    min_score = n_col ;
  1.2138 +    /* place in reverse order, so low column indices are at the front */
  1.2139 +    /* of the lists.  This is to encourage natural tie-breaking */
  1.2140 +    for (c = n_col-1 ; c >= 0 ; c--)
  1.2141 +    {
  1.2142 +        /* only add principal columns to degree lists */
  1.2143 +        if (COL_IS_ALIVE (c))
  1.2144 +        {
  1.2145 +            DEBUG4 (("place %d score %d minscore %d ncol %d\n",
  1.2146 +                c, Col [c].shared2.score, min_score, n_col)) ;
  1.2147 +
  1.2148 +            /* === Add columns score to DList =============================== */
  1.2149 +
  1.2150 +            score = Col [c].shared2.score ;
  1.2151 +
  1.2152 +            ASSERT (min_score >= 0) ;
  1.2153 +            ASSERT (min_score <= n_col) ;
  1.2154 +            ASSERT (score >= 0) ;
  1.2155 +            ASSERT (score <= n_col) ;
  1.2156 +            ASSERT (head [score] >= EMPTY) ;
  1.2157 +
  1.2158 +            /* now add this column to dList at proper score location */
  1.2159 +            next_col = head [score] ;
  1.2160 +            Col [c].shared3.prev = EMPTY ;
  1.2161 +            Col [c].shared4.degree_next = next_col ;
  1.2162 +
  1.2163 +            /* if there already was a column with the same score, set its */
  1.2164 +            /* previous pointer to this new column */
  1.2165 +            if (next_col != EMPTY)
  1.2166 +            {
  1.2167 +                Col [next_col].shared3.prev = c ;
  1.2168 +            }
  1.2169 +            head [score] = c ;
  1.2170 +
  1.2171 +            /* see if this score is less than current min */
  1.2172 +            min_score = MIN (min_score, score) ;
  1.2173 +
  1.2174 +#ifndef NDEBUG
  1.2175 +            debug_count++ ;
  1.2176 +#endif /* NDEBUG */
  1.2177 +
  1.2178 +        }
  1.2179 +    }
  1.2180 +
  1.2181 +#ifndef NDEBUG
  1.2182 +    DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
  1.2183 +        debug_count, n_col, n_col-debug_count)) ;
  1.2184 +    ASSERT (debug_count == n_col2) ;
  1.2185 +    debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
  1.2186 +#endif /* NDEBUG */
  1.2187 +
  1.2188 +    /* === Return number of remaining columns, and max row degree =========== */
  1.2189 +
  1.2190 +    *p_n_col2 = n_col2 ;
  1.2191 +    *p_n_row2 = n_row2 ;
  1.2192 +    *p_max_deg = max_deg ;
  1.2193 +}
  1.2194 +
  1.2195 +
  1.2196 +/* ========================================================================== */
  1.2197 +/* === find_ordering ======================================================== */
  1.2198 +/* ========================================================================== */
  1.2199 +
  1.2200 +/*
  1.2201 +    Order the principal columns of the supercolumn form of the matrix
  1.2202 +    (no supercolumns on input).  Uses a minimum approximate column minimum
  1.2203 +    degree ordering method.  Not user-callable.
  1.2204 +*/
  1.2205 +
  1.2206 +PRIVATE Int find_ordering       /* return the number of garbage collections */
  1.2207 +(
  1.2208 +    /* === Parameters ======================================================= */
  1.2209 +
  1.2210 +    Int n_row,                  /* number of rows of A */
  1.2211 +    Int n_col,                  /* number of columns of A */
  1.2212 +    Int Alen,                   /* size of A, 2*nnz + n_col or larger */
  1.2213 +    Colamd_Row Row [],          /* of size n_row+1 */
  1.2214 +    Colamd_Col Col [],          /* of size n_col+1 */
  1.2215 +    Int A [],                   /* column form and row form of A */
  1.2216 +    Int head [],                /* of size n_col+1 */
  1.2217 +    Int n_col2,                 /* Remaining columns to order */
  1.2218 +    Int max_deg,                /* Maximum row degree */
  1.2219 +    Int pfree,                  /* index of first free slot (2*nnz on entry) */
  1.2220 +    Int aggressive
  1.2221 +)
  1.2222 +{
  1.2223 +    /* === Local variables ================================================== */
  1.2224 +
  1.2225 +    Int k ;                     /* current pivot ordering step */
  1.2226 +    Int pivot_col ;             /* current pivot column */
  1.2227 +    Int *cp ;                   /* a column pointer */
  1.2228 +    Int *rp ;                   /* a row pointer */
  1.2229 +    Int pivot_row ;             /* current pivot row */
  1.2230 +    Int *new_cp ;               /* modified column pointer */
  1.2231 +    Int *new_rp ;               /* modified row pointer */
  1.2232 +    Int pivot_row_start ;       /* pointer to start of pivot row */
  1.2233 +    Int pivot_row_degree ;      /* number of columns in pivot row */
  1.2234 +    Int pivot_row_length ;      /* number of supercolumns in pivot row */
  1.2235 +    Int pivot_col_score ;       /* score of pivot column */
  1.2236 +    Int needed_memory ;         /* free space needed for pivot row */
  1.2237 +    Int *cp_end ;               /* pointer to the end of a column */
  1.2238 +    Int *rp_end ;               /* pointer to the end of a row */
  1.2239 +    Int row ;                   /* a row index */
  1.2240 +    Int col ;                   /* a column index */
  1.2241 +    Int max_score ;             /* maximum possible score */
  1.2242 +    Int cur_score ;             /* score of current column */
  1.2243 +    unsigned Int hash ;         /* hash value for supernode detection */
  1.2244 +    Int head_column ;           /* head of hash bucket */
  1.2245 +    Int first_col ;             /* first column in hash bucket */
  1.2246 +    Int tag_mark ;              /* marker value for mark array */
  1.2247 +    Int row_mark ;              /* Row [row].shared2.mark */
  1.2248 +    Int set_difference ;        /* set difference size of row with pivot row */
  1.2249 +    Int min_score ;             /* smallest column score */
  1.2250 +    Int col_thickness ;         /* "thickness" (no. of columns in a supercol) */
  1.2251 +    Int max_mark ;              /* maximum value of tag_mark */
  1.2252 +    Int pivot_col_thickness ;   /* number of columns represented by pivot col */
  1.2253 +    Int prev_col ;              /* Used by Dlist operations. */
  1.2254 +    Int next_col ;              /* Used by Dlist operations. */
  1.2255 +    Int ngarbage ;              /* number of garbage collections performed */
  1.2256 +
  1.2257 +#ifndef NDEBUG
  1.2258 +    Int debug_d ;               /* debug loop counter */
  1.2259 +    Int debug_step = 0 ;        /* debug loop counter */
  1.2260 +#endif /* NDEBUG */
  1.2261 +
  1.2262 +    /* === Initialization and clear mark ==================================== */
  1.2263 +
  1.2264 +    max_mark = INT_MAX - n_col ;        /* INT_MAX defined in <limits.h> */
  1.2265 +    tag_mark = clear_mark (0, max_mark, n_row, Row) ;
  1.2266 +    min_score = 0 ;
  1.2267 +    ngarbage = 0 ;
  1.2268 +    DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
  1.2269 +
  1.2270 +    /* === Order the columns ================================================ */
  1.2271 +
  1.2272 +    for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
  1.2273 +    {
  1.2274 +
  1.2275 +#ifndef NDEBUG
  1.2276 +        if (debug_step % 100 == 0)
  1.2277 +        {
  1.2278 +            DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
  1.2279 +        }
  1.2280 +        else
  1.2281 +        {
  1.2282 +            DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
  1.2283 +        }
  1.2284 +        debug_step++ ;
  1.2285 +        debug_deg_lists (n_row, n_col, Row, Col, head,
  1.2286 +                min_score, n_col2-k, max_deg) ;
  1.2287 +        debug_matrix (n_row, n_col, Row, Col, A) ;
  1.2288 +#endif /* NDEBUG */
  1.2289 +
  1.2290 +        /* === Select pivot column, and order it ============================ */
  1.2291 +
  1.2292 +        /* make sure degree list isn't empty */
  1.2293 +        ASSERT (min_score >= 0) ;
  1.2294 +        ASSERT (min_score <= n_col) ;
  1.2295 +        ASSERT (head [min_score] >= EMPTY) ;
  1.2296 +
  1.2297 +#ifndef NDEBUG
  1.2298 +        for (debug_d = 0 ; debug_d < min_score ; debug_d++)
  1.2299 +        {
  1.2300 +            ASSERT (head [debug_d] == EMPTY) ;
  1.2301 +        }
  1.2302 +#endif /* NDEBUG */
  1.2303 +
  1.2304 +        /* get pivot column from head of minimum degree list */
  1.2305 +        while (head [min_score] == EMPTY && min_score < n_col)
  1.2306 +        {
  1.2307 +            min_score++ ;
  1.2308 +        }
  1.2309 +        pivot_col = head [min_score] ;
  1.2310 +        ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
  1.2311 +        next_col = Col [pivot_col].shared4.degree_next ;
  1.2312 +        head [min_score] = next_col ;
  1.2313 +        if (next_col != EMPTY)
  1.2314 +        {
  1.2315 +            Col [next_col].shared3.prev = EMPTY ;
  1.2316 +        }
  1.2317 +
  1.2318 +        ASSERT (COL_IS_ALIVE (pivot_col)) ;
  1.2319 +
  1.2320 +        /* remember score for defrag check */
  1.2321 +        pivot_col_score = Col [pivot_col].shared2.score ;
  1.2322 +
  1.2323 +        /* the pivot column is the kth column in the pivot order */
  1.2324 +        Col [pivot_col].shared2.order = k ;
  1.2325 +
  1.2326 +        /* increment order count by column thickness */
  1.2327 +        pivot_col_thickness = Col [pivot_col].shared1.thickness ;
  1.2328 +        k += pivot_col_thickness ;
  1.2329 +        ASSERT (pivot_col_thickness > 0) ;
  1.2330 +        DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
  1.2331 +
  1.2332 +        /* === Garbage_collection, if necessary ============================= */
  1.2333 +
  1.2334 +        needed_memory = MIN (pivot_col_score, n_col - k) ;
  1.2335 +        if (pfree + needed_memory >= Alen)
  1.2336 +        {
  1.2337 +            pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
  1.2338 +            ngarbage++ ;
  1.2339 +            /* after garbage collection we will have enough */
  1.2340 +            ASSERT (pfree + needed_memory < Alen) ;
  1.2341 +            /* garbage collection has wiped out the Row[].shared2.mark array */
  1.2342 +            tag_mark = clear_mark (0, max_mark, n_row, Row) ;
  1.2343 +
  1.2344 +#ifndef NDEBUG
  1.2345 +            debug_matrix (n_row, n_col, Row, Col, A) ;
  1.2346 +#endif /* NDEBUG */
  1.2347 +        }
  1.2348 +
  1.2349 +        /* === Compute pivot row pattern ==================================== */
  1.2350 +
  1.2351 +        /* get starting location for this new merged row */
  1.2352 +        pivot_row_start = pfree ;
  1.2353 +
  1.2354 +        /* initialize new row counts to zero */
  1.2355 +        pivot_row_degree = 0 ;
  1.2356 +
  1.2357 +        /* tag pivot column as having been visited so it isn't included */
  1.2358 +        /* in merged pivot row */
  1.2359 +        Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
  1.2360 +
  1.2361 +        /* pivot row is the union of all rows in the pivot column pattern */
  1.2362 +        cp = &A [Col [pivot_col].start] ;
  1.2363 +        cp_end = cp + Col [pivot_col].length ;
  1.2364 +        while (cp < cp_end)
  1.2365 +        {
  1.2366 +            /* get a row */
  1.2367 +            row = *cp++ ;
  1.2368 +            DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
  1.2369 +            /* skip if row is dead */
  1.2370 +            if (ROW_IS_ALIVE (row))
  1.2371 +            {
  1.2372 +                rp = &A [Row [row].start] ;
  1.2373 +                rp_end = rp + Row [row].length ;
  1.2374 +                while (rp < rp_end)
  1.2375 +                {
  1.2376 +                    /* get a column */
  1.2377 +                    col = *rp++ ;
  1.2378 +                    /* add the column, if alive and untagged */
  1.2379 +                    col_thickness = Col [col].shared1.thickness ;
  1.2380 +                    if (col_thickness > 0 && COL_IS_ALIVE (col))
  1.2381 +                    {
  1.2382 +                        /* tag column in pivot row */
  1.2383 +                        Col [col].shared1.thickness = -col_thickness ;
  1.2384 +                        ASSERT (pfree < Alen) ;
  1.2385 +                        /* place column in pivot row */
  1.2386 +                        A [pfree++] = col ;
  1.2387 +                        pivot_row_degree += col_thickness ;
  1.2388 +                    }
  1.2389 +                }
  1.2390 +            }
  1.2391 +        }
  1.2392 +
  1.2393 +        /* clear tag on pivot column */
  1.2394 +        Col [pivot_col].shared1.thickness = pivot_col_thickness ;
  1.2395 +        max_deg = MAX (max_deg, pivot_row_degree) ;
  1.2396 +
  1.2397 +#ifndef NDEBUG
  1.2398 +        DEBUG3 (("check2\n")) ;
  1.2399 +        debug_mark (n_row, Row, tag_mark, max_mark) ;
  1.2400 +#endif /* NDEBUG */
  1.2401 +
  1.2402 +        /* === Kill all rows used to construct pivot row ==================== */
  1.2403 +
  1.2404 +        /* also kill pivot row, temporarily */
  1.2405 +        cp = &A [Col [pivot_col].start] ;
  1.2406 +        cp_end = cp + Col [pivot_col].length ;
  1.2407 +        while (cp < cp_end)
  1.2408 +        {
  1.2409 +            /* may be killing an already dead row */
  1.2410 +            row = *cp++ ;
  1.2411 +            DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
  1.2412 +            KILL_ROW (row) ;
  1.2413 +        }
  1.2414 +
  1.2415 +        /* === Select a row index to use as the new pivot row =============== */
  1.2416 +
  1.2417 +        pivot_row_length = pfree - pivot_row_start ;
  1.2418 +        if (pivot_row_length > 0)
  1.2419 +        {
  1.2420 +            /* pick the "pivot" row arbitrarily (first row in col) */
  1.2421 +            pivot_row = A [Col [pivot_col].start] ;
  1.2422 +            DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
  1.2423 +        }
  1.2424 +        else
  1.2425 +        {
  1.2426 +            /* there is no pivot row, since it is of zero length */
  1.2427 +            pivot_row = EMPTY ;
  1.2428 +            ASSERT (pivot_row_length == 0) ;
  1.2429 +        }
  1.2430 +        ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
  1.2431 +
  1.2432 +        /* === Approximate degree computation =============================== */
  1.2433 +
  1.2434 +        /* Here begins the computation of the approximate degree.  The column */
  1.2435 +        /* score is the sum of the pivot row "length", plus the size of the */
  1.2436 +        /* set differences of each row in the column minus the pattern of the */
  1.2437 +        /* pivot row itself.  The column ("thickness") itself is also */
  1.2438 +        /* excluded from the column score (we thus use an approximate */
  1.2439 +        /* external degree). */
  1.2440 +
  1.2441 +        /* The time taken by the following code (compute set differences, and */
  1.2442 +        /* add them up) is proportional to the size of the data structure */
  1.2443 +        /* being scanned - that is, the sum of the sizes of each column in */
  1.2444 +        /* the pivot row.  Thus, the amortized time to compute a column score */
  1.2445 +        /* is proportional to the size of that column (where size, in this */
  1.2446 +        /* context, is the column "length", or the number of row indices */
  1.2447 +        /* in that column).  The number of row indices in a column is */
  1.2448 +        /* monotonically non-decreasing, from the length of the original */
  1.2449 +        /* column on input to colamd. */
  1.2450 +
  1.2451 +        /* === Compute set differences ====================================== */
  1.2452 +
  1.2453 +        DEBUG3 (("** Computing set differences phase. **\n")) ;
  1.2454 +
  1.2455 +        /* pivot row is currently dead - it will be revived later. */
  1.2456 +
  1.2457 +        DEBUG3 (("Pivot row: ")) ;
  1.2458 +        /* for each column in pivot row */
  1.2459 +        rp = &A [pivot_row_start] ;
  1.2460 +        rp_end = rp + pivot_row_length ;
  1.2461 +        while (rp < rp_end)
  1.2462 +        {
  1.2463 +            col = *rp++ ;
  1.2464 +            ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
  1.2465 +            DEBUG3 (("Col: %d\n", col)) ;
  1.2466 +
  1.2467 +            /* clear tags used to construct pivot row pattern */
  1.2468 +            col_thickness = -Col [col].shared1.thickness ;
  1.2469 +            ASSERT (col_thickness > 0) ;
  1.2470 +            Col [col].shared1.thickness = col_thickness ;
  1.2471 +
  1.2472 +            /* === Remove column from degree list =========================== */
  1.2473 +
  1.2474 +            cur_score = Col [col].shared2.score ;
  1.2475 +            prev_col = Col [col].shared3.prev ;
  1.2476 +            next_col = Col [col].shared4.degree_next ;
  1.2477 +            ASSERT (cur_score >= 0) ;
  1.2478 +            ASSERT (cur_score <= n_col) ;
  1.2479 +            ASSERT (cur_score >= EMPTY) ;
  1.2480 +            if (prev_col == EMPTY)
  1.2481 +            {
  1.2482 +                head [cur_score] = next_col ;
  1.2483 +            }
  1.2484 +            else
  1.2485 +            {
  1.2486 +                Col [prev_col].shared4.degree_next = next_col ;
  1.2487 +            }
  1.2488 +            if (next_col != EMPTY)
  1.2489 +            {
  1.2490 +                Col [next_col].shared3.prev = prev_col ;
  1.2491 +            }
  1.2492 +
  1.2493 +            /* === Scan the column ========================================== */
  1.2494 +
  1.2495 +            cp = &A [Col [col].start] ;
  1.2496 +            cp_end = cp + Col [col].length ;
  1.2497 +            while (cp < cp_end)
  1.2498 +            {
  1.2499 +                /* get a row */
  1.2500 +                row = *cp++ ;
  1.2501 +                row_mark = Row [row].shared2.mark ;
  1.2502 +                /* skip if dead */
  1.2503 +                if (ROW_IS_MARKED_DEAD (row_mark))
  1.2504 +                {
  1.2505 +                    continue ;
  1.2506 +                }
  1.2507 +                ASSERT (row != pivot_row) ;
  1.2508 +                set_difference = row_mark - tag_mark ;
  1.2509 +                /* check if the row has been seen yet */
  1.2510 +                if (set_difference < 0)
  1.2511 +                {
  1.2512 +                    ASSERT (Row [row].shared1.degree <= max_deg) ;
  1.2513 +                    set_difference = Row [row].shared1.degree ;
  1.2514 +                }
  1.2515 +                /* subtract column thickness from this row's set difference */
  1.2516 +                set_difference -= col_thickness ;
  1.2517 +                ASSERT (set_difference >= 0) ;
  1.2518 +                /* absorb this row if the set difference becomes zero */
  1.2519 +                if (set_difference == 0 && aggressive)
  1.2520 +                {
  1.2521 +                    DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
  1.2522 +                    KILL_ROW (row) ;
  1.2523 +                }
  1.2524 +                else
  1.2525 +                {
  1.2526 +                    /* save the new mark */
  1.2527 +                    Row [row].shared2.mark = set_difference + tag_mark ;
  1.2528 +                }
  1.2529 +            }
  1.2530 +        }
  1.2531 +
  1.2532 +#ifndef NDEBUG
  1.2533 +        debug_deg_lists (n_row, n_col, Row, Col, head,
  1.2534 +                min_score, n_col2-k-pivot_row_degree, max_deg) ;
  1.2535 +#endif /* NDEBUG */
  1.2536 +
  1.2537 +        /* === Add up set differences for each column ======================= */
  1.2538 +
  1.2539 +        DEBUG3 (("** Adding set differences phase. **\n")) ;
  1.2540 +
  1.2541 +        /* for each column in pivot row */
  1.2542 +        rp = &A [pivot_row_start] ;
  1.2543 +        rp_end = rp + pivot_row_length ;
  1.2544 +        while (rp < rp_end)
  1.2545 +        {
  1.2546 +            /* get a column */
  1.2547 +            col = *rp++ ;
  1.2548 +            ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
  1.2549 +            hash = 0 ;
  1.2550 +            cur_score = 0 ;
  1.2551 +            cp = &A [Col [col].start] ;
  1.2552 +            /* compact the column */
  1.2553 +            new_cp = cp ;
  1.2554 +            cp_end = cp + Col [col].length ;
  1.2555 +
  1.2556 +            DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
  1.2557 +
  1.2558 +            while (cp < cp_end)
  1.2559 +            {
  1.2560 +                /* get a row */
  1.2561 +                row = *cp++ ;
  1.2562 +                ASSERT(row >= 0 && row < n_row) ;
  1.2563 +                row_mark = Row [row].shared2.mark ;
  1.2564 +                /* skip if dead */
  1.2565 +                if (ROW_IS_MARKED_DEAD (row_mark))
  1.2566 +                {
  1.2567 +                    DEBUG4 ((" Row %d, dead\n", row)) ;
  1.2568 +                    continue ;
  1.2569 +                }
  1.2570 +                DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
  1.2571 +                ASSERT (row_mark >= tag_mark) ;
  1.2572 +                /* compact the column */
  1.2573 +                *new_cp++ = row ;
  1.2574 +                /* compute hash function */
  1.2575 +                hash += row ;
  1.2576 +                /* add set difference */
  1.2577 +                cur_score += row_mark - tag_mark ;
  1.2578 +                /* integer overflow... */
  1.2579 +                cur_score = MIN (cur_score, n_col) ;
  1.2580 +            }
  1.2581 +
  1.2582 +            /* recompute the column's length */
  1.2583 +            Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
  1.2584 +
  1.2585 +            /* === Further mass elimination ================================= */
  1.2586 +
  1.2587 +            if (Col [col].length == 0)
  1.2588 +            {
  1.2589 +                DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
  1.2590 +                /* nothing left but the pivot row in this column */
  1.2591 +                KILL_PRINCIPAL_COL (col) ;
  1.2592 +                pivot_row_degree -= Col [col].shared1.thickness ;
  1.2593 +                ASSERT (pivot_row_degree >= 0) ;
  1.2594 +                /* order it */
  1.2595 +                Col [col].shared2.order = k ;
  1.2596 +                /* increment order count by column thickness */
  1.2597 +                k += Col [col].shared1.thickness ;
  1.2598 +            }
  1.2599 +            else
  1.2600 +            {
  1.2601 +                /* === Prepare for supercolumn detection ==================== */
  1.2602 +
  1.2603 +                DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
  1.2604 +
  1.2605 +                /* save score so far */
  1.2606 +                Col [col].shared2.score = cur_score ;
  1.2607 +
  1.2608 +                /* add column to hash table, for supercolumn detection */
  1.2609 +                hash %= n_col + 1 ;
  1.2610 +
  1.2611 +                DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
  1.2612 +                ASSERT (((Int) hash) <= n_col) ;
  1.2613 +
  1.2614 +                head_column = head [hash] ;
  1.2615 +                if (head_column > EMPTY)
  1.2616 +                {
  1.2617 +                    /* degree list "hash" is non-empty, use prev (shared3) of */
  1.2618 +                    /* first column in degree list as head of hash bucket */
  1.2619 +                    first_col = Col [head_column].shared3.headhash ;
  1.2620 +                    Col [head_column].shared3.headhash = col ;
  1.2621 +                }
  1.2622 +                else
  1.2623 +                {
  1.2624 +                    /* degree list "hash" is empty, use head as hash bucket */
  1.2625 +                    first_col = - (head_column + 2) ;
  1.2626 +                    head [hash] = - (col + 2) ;
  1.2627 +                }
  1.2628 +                Col [col].shared4.hash_next = first_col ;
  1.2629 +
  1.2630 +                /* save hash function in Col [col].shared3.hash */
  1.2631 +                Col [col].shared3.hash = (Int) hash ;
  1.2632 +                ASSERT (COL_IS_ALIVE (col)) ;
  1.2633 +            }
  1.2634 +        }
  1.2635 +
  1.2636 +        /* The approximate external column degree is now computed.  */
  1.2637 +
  1.2638 +        /* === Supercolumn detection ======================================== */
  1.2639 +
  1.2640 +        DEBUG3 (("** Supercolumn detection phase. **\n")) ;
  1.2641 +
  1.2642 +        detect_super_cols (
  1.2643 +
  1.2644 +#ifndef NDEBUG
  1.2645 +                n_col, Row,
  1.2646 +#endif /* NDEBUG */
  1.2647 +
  1.2648 +                Col, A, head, pivot_row_start, pivot_row_length) ;
  1.2649 +
  1.2650 +        /* === Kill the pivotal column ====================================== */
  1.2651 +
  1.2652 +        KILL_PRINCIPAL_COL (pivot_col) ;
  1.2653 +
  1.2654 +        /* === Clear mark =================================================== */
  1.2655 +
  1.2656 +        tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
  1.2657 +
  1.2658 +#ifndef NDEBUG
  1.2659 +        DEBUG3 (("check3\n")) ;
  1.2660 +        debug_mark (n_row, Row, tag_mark, max_mark) ;
  1.2661 +#endif /* NDEBUG */
  1.2662 +
  1.2663 +        /* === Finalize the new pivot row, and column scores ================ */
  1.2664 +
  1.2665 +        DEBUG3 (("** Finalize scores phase. **\n")) ;
  1.2666 +
  1.2667 +        /* for each column in pivot row */
  1.2668 +        rp = &A [pivot_row_start] ;
  1.2669 +        /* compact the pivot row */
  1.2670 +        new_rp = rp ;
  1.2671 +        rp_end = rp + pivot_row_length ;
  1.2672 +        while (rp < rp_end)
  1.2673 +        {
  1.2674 +            col = *rp++ ;
  1.2675 +            /* skip dead columns */
  1.2676 +            if (COL_IS_DEAD (col))
  1.2677 +            {
  1.2678 +                continue ;
  1.2679 +            }
  1.2680 +            *new_rp++ = col ;
  1.2681 +            /* add new pivot row to column */
  1.2682 +            A [Col [col].start + (Col [col].length++)] = pivot_row ;
  1.2683 +
  1.2684 +            /* retrieve score so far and add on pivot row's degree. */
  1.2685 +            /* (we wait until here for this in case the pivot */
  1.2686 +            /* row's degree was reduced due to mass elimination). */
  1.2687 +            cur_score = Col [col].shared2.score + pivot_row_degree ;
  1.2688 +
  1.2689 +            /* calculate the max possible score as the number of */
  1.2690 +            /* external columns minus the 'k' value minus the */
  1.2691 +            /* columns thickness */
  1.2692 +            max_score = n_col - k - Col [col].shared1.thickness ;
  1.2693 +
  1.2694 +            /* make the score the external degree of the union-of-rows */
  1.2695 +            cur_score -= Col [col].shared1.thickness ;
  1.2696 +
  1.2697 +            /* make sure score is less or equal than the max score */
  1.2698 +            cur_score = MIN (cur_score, max_score) ;
  1.2699 +            ASSERT (cur_score >= 0) ;
  1.2700 +
  1.2701 +            /* store updated score */
  1.2702 +            Col [col].shared2.score = cur_score ;
  1.2703 +
  1.2704 +            /* === Place column back in degree list ========================= */
  1.2705 +
  1.2706 +            ASSERT (min_score >= 0) ;
  1.2707 +            ASSERT (min_score <= n_col) ;
  1.2708 +            ASSERT (cur_score >= 0) ;
  1.2709 +            ASSERT (cur_score <= n_col) ;
  1.2710 +            ASSERT (head [cur_score] >= EMPTY) ;
  1.2711 +            next_col = head [cur_score] ;
  1.2712 +            Col [col].shared4.degree_next = next_col ;
  1.2713 +            Col [col].shared3.prev = EMPTY ;
  1.2714 +            if (next_col != EMPTY)
  1.2715 +            {
  1.2716 +                Col [next_col].shared3.prev = col ;
  1.2717 +            }
  1.2718 +            head [cur_score] = col ;
  1.2719 +
  1.2720 +            /* see if this score is less than current min */
  1.2721 +            min_score = MIN (min_score, cur_score) ;
  1.2722 +
  1.2723 +        }
  1.2724 +
  1.2725 +#ifndef NDEBUG
  1.2726 +        debug_deg_lists (n_row, n_col, Row, Col, head,
  1.2727 +                min_score, n_col2-k, max_deg) ;
  1.2728 +#endif /* NDEBUG */
  1.2729 +
  1.2730 +        /* === Resurrect the new pivot row ================================== */
  1.2731 +
  1.2732 +        if (pivot_row_degree > 0)
  1.2733 +        {
  1.2734 +            /* update pivot row length to reflect any cols that were killed */
  1.2735 +            /* during super-col detection and mass elimination */
  1.2736 +            Row [pivot_row].start  = pivot_row_start ;
  1.2737 +            Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
  1.2738 +            ASSERT (Row [pivot_row].length > 0) ;
  1.2739 +            Row [pivot_row].shared1.degree = pivot_row_degree ;
  1.2740 +            Row [pivot_row].shared2.mark = 0 ;
  1.2741 +            /* pivot row is no longer dead */
  1.2742 +
  1.2743 +            DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
  1.2744 +                        pivot_row, pivot_row_degree)) ;
  1.2745 +        }
  1.2746 +    }
  1.2747 +
  1.2748 +    /* === All principal columns have now been ordered ====================== */
  1.2749 +
  1.2750 +    return (ngarbage) ;
  1.2751 +}
  1.2752 +
  1.2753 +
  1.2754 +/* ========================================================================== */
  1.2755 +/* === order_children ======================================================= */
  1.2756 +/* ========================================================================== */
  1.2757 +
  1.2758 +/*
  1.2759 +    The find_ordering routine has ordered all of the principal columns (the
  1.2760 +    representatives of the supercolumns).  The non-principal columns have not
  1.2761 +    yet been ordered.  This routine orders those columns by walking up the
  1.2762 +    parent tree (a column is a child of the column which absorbed it).  The
  1.2763 +    final permutation vector is then placed in p [0 ... n_col-1], with p [0]
  1.2764 +    being the first column, and p [n_col-1] being the last.  It doesn't look
  1.2765 +    like it at first glance, but be assured that this routine takes time linear
  1.2766 +    in the number of columns.  Although not immediately obvious, the time
  1.2767 +    taken by this routine is O (n_col), that is, linear in the number of
  1.2768 +    columns.  Not user-callable.
  1.2769 +*/
  1.2770 +
  1.2771 +PRIVATE void order_children
  1.2772 +(
  1.2773 +    /* === Parameters ======================================================= */
  1.2774 +
  1.2775 +    Int n_col,                  /* number of columns of A */
  1.2776 +    Colamd_Col Col [],          /* of size n_col+1 */
  1.2777 +    Int p []                    /* p [0 ... n_col-1] is the column permutation*/
  1.2778 +)
  1.2779 +{
  1.2780 +    /* === Local variables ================================================== */
  1.2781 +
  1.2782 +    Int i ;                     /* loop counter for all columns */
  1.2783 +    Int c ;                     /* column index */
  1.2784 +    Int parent ;                /* index of column's parent */
  1.2785 +    Int order ;                 /* column's order */
  1.2786 +
  1.2787 +    /* === Order each non-principal column ================================== */
  1.2788 +
  1.2789 +    for (i = 0 ; i < n_col ; i++)
  1.2790 +    {
  1.2791 +        /* find an un-ordered non-principal column */
  1.2792 +        ASSERT (COL_IS_DEAD (i)) ;
  1.2793 +        if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
  1.2794 +        {
  1.2795 +            parent = i ;
  1.2796 +            /* once found, find its principal parent */
  1.2797 +            do
  1.2798 +            {
  1.2799 +                parent = Col [parent].shared1.parent ;
  1.2800 +            } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
  1.2801 +
  1.2802 +            /* now, order all un-ordered non-principal columns along path */
  1.2803 +            /* to this parent.  collapse tree at the same time */
  1.2804 +            c = i ;
  1.2805 +            /* get order of parent */
  1.2806 +            order = Col [parent].shared2.order ;
  1.2807 +
  1.2808 +            do
  1.2809 +            {
  1.2810 +                ASSERT (Col [c].shared2.order == EMPTY) ;
  1.2811 +
  1.2812 +                /* order this column */
  1.2813 +                Col [c].shared2.order = order++ ;
  1.2814 +                /* collaps tree */
  1.2815 +                Col [c].shared1.parent = parent ;
  1.2816 +
  1.2817 +                /* get immediate parent of this column */
  1.2818 +                c = Col [c].shared1.parent ;
  1.2819 +
  1.2820 +                /* continue until we hit an ordered column.  There are */
  1.2821 +                /* guarranteed not to be anymore unordered columns */
  1.2822 +                /* above an ordered column */
  1.2823 +            } while (Col [c].shared2.order == EMPTY) ;
  1.2824 +
  1.2825 +            /* re-order the super_col parent to largest order for this group */
  1.2826 +            Col [parent].shared2.order = order ;
  1.2827 +        }
  1.2828 +    }
  1.2829 +
  1.2830 +    /* === Generate the permutation ========================================= */
  1.2831 +
  1.2832 +    for (c = 0 ; c < n_col ; c++)
  1.2833 +    {
  1.2834 +        p [Col [c].shared2.order] = c ;
  1.2835 +    }
  1.2836 +}
  1.2837 +
  1.2838 +
  1.2839 +/* ========================================================================== */
  1.2840 +/* === detect_super_cols ==================================================== */
  1.2841 +/* ========================================================================== */
  1.2842 +
  1.2843 +/*
  1.2844 +    Detects supercolumns by finding matches between columns in the hash buckets.
  1.2845 +    Check amongst columns in the set A [row_start ... row_start + row_length-1].
  1.2846 +    The columns under consideration are currently *not* in the degree lists,
  1.2847 +    and have already been placed in the hash buckets.
  1.2848 +
  1.2849 +    The hash bucket for columns whose hash function is equal to h is stored
  1.2850 +    as follows:
  1.2851 +
  1.2852 +        if head [h] is >= 0, then head [h] contains a degree list, so:
  1.2853 +
  1.2854 +                head [h] is the first column in degree bucket h.
  1.2855 +                Col [head [h]].headhash gives the first column in hash bucket h.
  1.2856 +
  1.2857 +        otherwise, the degree list is empty, and:
  1.2858 +
  1.2859 +                -(head [h] + 2) is the first column in hash bucket h.
  1.2860 +
  1.2861 +    For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
  1.2862 +    column" pointer.  Col [c].shared3.hash is used instead as the hash number
  1.2863 +    for that column.  The value of Col [c].shared4.hash_next is the next column
  1.2864 +    in the same hash bucket.
  1.2865 +
  1.2866 +    Assuming no, or "few" hash collisions, the time taken by this routine is
  1.2867 +    linear in the sum of the sizes (lengths) of each column whose score has
  1.2868 +    just been computed in the approximate degree computation.
  1.2869 +    Not user-callable.
  1.2870 +*/
  1.2871 +
  1.2872 +PRIVATE void detect_super_cols
  1.2873 +(
  1.2874 +    /* === Parameters ======================================================= */
  1.2875 +
  1.2876 +#ifndef NDEBUG
  1.2877 +    /* these two parameters are only needed when debugging is enabled: */
  1.2878 +    Int n_col,                  /* number of columns of A */
  1.2879 +    Colamd_Row Row [],          /* of size n_row+1 */
  1.2880 +#endif /* NDEBUG */
  1.2881 +
  1.2882 +    Colamd_Col Col [],          /* of size n_col+1 */
  1.2883 +    Int A [],                   /* row indices of A */
  1.2884 +    Int head [],                /* head of degree lists and hash buckets */
  1.2885 +    Int row_start,              /* pointer to set of columns to check */
  1.2886 +    Int row_length              /* number of columns to check */
  1.2887 +)
  1.2888 +{
  1.2889 +    /* === Local variables ================================================== */
  1.2890 +
  1.2891 +    Int hash ;                  /* hash value for a column */
  1.2892 +    Int *rp ;                   /* pointer to a row */
  1.2893 +    Int c ;                     /* a column index */
  1.2894 +    Int super_c ;               /* column index of the column to absorb into */
  1.2895 +    Int *cp1 ;                  /* column pointer for column super_c */
  1.2896 +    Int *cp2 ;                  /* column pointer for column c */
  1.2897 +    Int length ;                /* length of column super_c */
  1.2898 +    Int prev_c ;                /* column preceding c in hash bucket */
  1.2899 +    Int i ;                     /* loop counter */
  1.2900 +    Int *rp_end ;               /* pointer to the end of the row */
  1.2901 +    Int col ;                   /* a column index in the row to check */
  1.2902 +    Int head_column ;           /* first column in hash bucket or degree list */
  1.2903 +    Int first_col ;             /* first column in hash bucket */
  1.2904 +
  1.2905 +    /* === Consider each column in the row ================================== */
  1.2906 +
  1.2907 +    rp = &A [row_start] ;
  1.2908 +    rp_end = rp + row_length ;
  1.2909 +    while (rp < rp_end)
  1.2910 +    {
  1.2911 +        col = *rp++ ;
  1.2912 +        if (COL_IS_DEAD (col))
  1.2913 +        {
  1.2914 +            continue ;
  1.2915 +        }
  1.2916 +
  1.2917 +        /* get hash number for this column */
  1.2918 +        hash = Col [col].shared3.hash ;
  1.2919 +        ASSERT (hash <= n_col) ;
  1.2920 +
  1.2921 +        /* === Get the first column in this hash bucket ===================== */
  1.2922 +
  1.2923 +        head_column = head [hash] ;
  1.2924 +        if (head_column > EMPTY)
  1.2925 +        {
  1.2926 +            first_col = Col [head_column].shared3.headhash ;
  1.2927 +        }
  1.2928 +        else
  1.2929 +        {
  1.2930 +            first_col = - (head_column + 2) ;
  1.2931 +        }
  1.2932 +
  1.2933 +        /* === Consider each column in the hash bucket ====================== */
  1.2934 +
  1.2935 +        for (super_c = first_col ; super_c != EMPTY ;
  1.2936 +            super_c = Col [super_c].shared4.hash_next)
  1.2937 +        {
  1.2938 +            ASSERT (COL_IS_ALIVE (super_c)) ;
  1.2939 +            ASSERT (Col [super_c].shared3.hash == hash) ;
  1.2940 +            length = Col [super_c].length ;
  1.2941 +
  1.2942 +            /* prev_c is the column preceding column c in the hash bucket */
  1.2943 +            prev_c = super_c ;
  1.2944 +
  1.2945 +            /* === Compare super_c with all columns after it ================ */
  1.2946 +
  1.2947 +            for (c = Col [super_c].shared4.hash_next ;
  1.2948 +                 c != EMPTY ; c = Col [c].shared4.hash_next)
  1.2949 +            {
  1.2950 +                ASSERT (c != super_c) ;
  1.2951 +                ASSERT (COL_IS_ALIVE (c)) ;
  1.2952 +                ASSERT (Col [c].shared3.hash == hash) ;
  1.2953 +
  1.2954 +                /* not identical if lengths or scores are different */
  1.2955 +                if (Col [c].length != length ||
  1.2956 +                    Col [c].shared2.score != Col [super_c].shared2.score)
  1.2957 +                {
  1.2958 +                    prev_c = c ;
  1.2959 +                    continue ;
  1.2960 +                }
  1.2961 +
  1.2962 +                /* compare the two columns */
  1.2963 +                cp1 = &A [Col [super_c].start] ;
  1.2964 +                cp2 = &A [Col [c].start] ;
  1.2965 +
  1.2966 +                for (i = 0 ; i < length ; i++)
  1.2967 +                {
  1.2968 +                    /* the columns are "clean" (no dead rows) */
  1.2969 +                    ASSERT (ROW_IS_ALIVE (*cp1))  ;
  1.2970 +                    ASSERT (ROW_IS_ALIVE (*cp2))  ;
  1.2971 +                    /* row indices will same order for both supercols, */
  1.2972 +                    /* no gather scatter nessasary */
  1.2973 +                    if (*cp1++ != *cp2++)
  1.2974 +                    {
  1.2975 +                        break ;
  1.2976 +                    }
  1.2977 +                }
  1.2978 +
  1.2979 +                /* the two columns are different if the for-loop "broke" */
  1.2980 +                if (i != length)
  1.2981 +                {
  1.2982 +                    prev_c = c ;
  1.2983 +                    continue ;
  1.2984 +                }
  1.2985 +
  1.2986 +                /* === Got it!  two columns are identical =================== */
  1.2987 +
  1.2988 +                ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
  1.2989 +
  1.2990 +                Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
  1.2991 +                Col [c].shared1.parent = super_c ;
  1.2992 +                KILL_NON_PRINCIPAL_COL (c) ;
  1.2993 +                /* order c later, in order_children() */
  1.2994 +                Col [c].shared2.order = EMPTY ;
  1.2995 +                /* remove c from hash bucket */
  1.2996 +                Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
  1.2997 +            }
  1.2998 +        }
  1.2999 +
  1.3000 +        /* === Empty this hash bucket ======================================= */
  1.3001 +
  1.3002 +        if (head_column > EMPTY)
  1.3003 +        {
  1.3004 +            /* corresponding degree list "hash" is not empty */
  1.3005 +            Col [head_column].shared3.headhash = EMPTY ;
  1.3006 +        }
  1.3007 +        else
  1.3008 +        {
  1.3009 +            /* corresponding degree list "hash" is empty */
  1.3010 +            head [hash] = EMPTY ;
  1.3011 +        }
  1.3012 +    }
  1.3013 +}
  1.3014 +
  1.3015 +
  1.3016 +/* ========================================================================== */
  1.3017 +/* === garbage_collection =================================================== */
  1.3018 +/* ========================================================================== */
  1.3019 +
  1.3020 +/*
  1.3021 +    Defragments and compacts columns and rows in the workspace A.  Used when
  1.3022 +    all avaliable memory has been used while performing row merging.  Returns
  1.3023 +    the index of the first free position in A, after garbage collection.  The
  1.3024 +    time taken by this routine is linear is the size of the array A, which is
  1.3025 +    itself linear in the number of nonzeros in the input matrix.
  1.3026 +    Not user-callable.
  1.3027 +*/
  1.3028 +
  1.3029 +PRIVATE Int garbage_collection  /* returns the new value of pfree */
  1.3030 +(
  1.3031 +    /* === Parameters ======================================================= */
  1.3032 +
  1.3033 +    Int n_row,                  /* number of rows */
  1.3034 +    Int n_col,                  /* number of columns */
  1.3035 +    Colamd_Row Row [],          /* row info */
  1.3036 +    Colamd_Col Col [],          /* column info */
  1.3037 +    Int A [],                   /* A [0 ... Alen-1] holds the matrix */
  1.3038 +    Int *pfree                  /* &A [0] ... pfree is in use */
  1.3039 +)
  1.3040 +{
  1.3041 +    /* === Local variables ================================================== */
  1.3042 +
  1.3043 +    Int *psrc ;                 /* source pointer */
  1.3044 +    Int *pdest ;                /* destination pointer */
  1.3045 +    Int j ;                     /* counter */
  1.3046 +    Int r ;                     /* a row index */
  1.3047 +    Int c ;                     /* a column index */
  1.3048 +    Int length ;                /* length of a row or column */
  1.3049 +
  1.3050 +#ifndef NDEBUG
  1.3051 +    Int debug_rows ;
  1.3052 +    DEBUG2 (("Defrag..\n")) ;
  1.3053 +    for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
  1.3054 +    debug_rows = 0 ;
  1.3055 +#endif /* NDEBUG */
  1.3056 +
  1.3057 +    /* === Defragment the columns =========================================== */
  1.3058 +
  1.3059 +    pdest = &A[0] ;
  1.3060 +    for (c = 0 ; c < n_col ; c++)
  1.3061 +    {
  1.3062 +        if (COL_IS_ALIVE (c))
  1.3063 +        {
  1.3064 +            psrc = &A [Col [c].start] ;
  1.3065 +
  1.3066 +            /* move and compact the column */
  1.3067 +            ASSERT (pdest <= psrc) ;
  1.3068 +            Col [c].start = (Int) (pdest - &A [0]) ;
  1.3069 +            length = Col [c].length ;
  1.3070 +            for (j = 0 ; j < length ; j++)
  1.3071 +            {
  1.3072 +                r = *psrc++ ;
  1.3073 +                if (ROW_IS_ALIVE (r))
  1.3074 +                {
  1.3075 +                    *pdest++ = r ;
  1.3076 +                }
  1.3077 +            }
  1.3078 +            Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
  1.3079 +        }
  1.3080 +    }
  1.3081 +
  1.3082 +    /* === Prepare to defragment the rows =================================== */
  1.3083 +
  1.3084 +    for (r = 0 ; r < n_row ; r++)
  1.3085 +    {
  1.3086 +        if (ROW_IS_DEAD (r) || (Row [r].length == 0))
  1.3087 +        {
  1.3088 +            /* This row is already dead, or is of zero length.  Cannot compact
  1.3089 +             * a row of zero length, so kill it.  NOTE: in the current version,
  1.3090 +             * there are no zero-length live rows.  Kill the row (for the first
  1.3091 +             * time, or again) just to be safe. */
  1.3092 +            KILL_ROW (r) ;
  1.3093 +        }
  1.3094 +        else
  1.3095 +        {
  1.3096 +            /* save first column index in Row [r].shared2.first_column */
  1.3097 +            psrc = &A [Row [r].start] ;
  1.3098 +            Row [r].shared2.first_column = *psrc ;
  1.3099 +            ASSERT (ROW_IS_ALIVE (r)) ;
  1.3100 +            /* flag the start of the row with the one's complement of row */
  1.3101 +            *psrc = ONES_COMPLEMENT (r) ;
  1.3102 +#ifndef NDEBUG
  1.3103 +            debug_rows++ ;
  1.3104 +#endif /* NDEBUG */
  1.3105 +        }
  1.3106 +    }
  1.3107 +
  1.3108 +    /* === Defragment the rows ============================================== */
  1.3109 +
  1.3110 +    psrc = pdest ;
  1.3111 +    while (psrc < pfree)
  1.3112 +    {
  1.3113 +        /* find a negative number ... the start of a row */
  1.3114 +        if (*psrc++ < 0)
  1.3115 +        {
  1.3116 +            psrc-- ;
  1.3117 +            /* get the row index */
  1.3118 +            r = ONES_COMPLEMENT (*psrc) ;
  1.3119 +            ASSERT (r >= 0 && r < n_row) ;
  1.3120 +            /* restore first column index */
  1.3121 +            *psrc = Row [r].shared2.first_column ;
  1.3122 +            ASSERT (ROW_IS_ALIVE (r)) ;
  1.3123 +            ASSERT (Row [r].length > 0) ;
  1.3124 +            /* move and compact the row */
  1.3125 +            ASSERT (pdest <= psrc) ;
  1.3126 +            Row [r].start = (Int) (pdest - &A [0]) ;
  1.3127 +            length = Row [r].length ;
  1.3128 +            for (j = 0 ; j < length ; j++)
  1.3129 +            {
  1.3130 +                c = *psrc++ ;
  1.3131 +                if (COL_IS_ALIVE (c))
  1.3132 +                {
  1.3133 +                    *pdest++ = c ;
  1.3134 +                }
  1.3135 +            }
  1.3136 +            Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
  1.3137 +            ASSERT (Row [r].length > 0) ;
  1.3138 +#ifndef NDEBUG
  1.3139 +            debug_rows-- ;
  1.3140 +#endif /* NDEBUG */
  1.3141 +        }
  1.3142 +    }
  1.3143 +    /* ensure we found all the rows */
  1.3144 +    ASSERT (debug_rows == 0) ;
  1.3145 +
  1.3146 +    /* === Return the new value of pfree ==================================== */
  1.3147 +
  1.3148 +    return ((Int) (pdest - &A [0])) ;
  1.3149 +}
  1.3150 +
  1.3151 +
  1.3152 +/* ========================================================================== */
  1.3153 +/* === clear_mark =========================================================== */
  1.3154 +/* ========================================================================== */
  1.3155 +
  1.3156 +/*
  1.3157 +    Clears the Row [].shared2.mark array, and returns the new tag_mark.
  1.3158 +    Return value is the new tag_mark.  Not user-callable.
  1.3159 +*/
  1.3160 +
  1.3161 +PRIVATE Int clear_mark  /* return the new value for tag_mark */
  1.3162 +(
  1.3163 +    /* === Parameters ======================================================= */
  1.3164 +
  1.3165 +    Int tag_mark,       /* new value of tag_mark */
  1.3166 +    Int max_mark,       /* max allowed value of tag_mark */
  1.3167 +
  1.3168 +    Int n_row,          /* number of rows in A */
  1.3169 +    Colamd_Row Row []   /* Row [0 ... n_row-1].shared2.mark is set to zero */
  1.3170 +)
  1.3171 +{
  1.3172 +    /* === Local variables ================================================== */
  1.3173 +
  1.3174 +    Int r ;
  1.3175 +
  1.3176 +    if (tag_mark <= 0 || tag_mark >= max_mark)
  1.3177 +    {
  1.3178 +        for (r = 0 ; r < n_row ; r++)
  1.3179 +        {
  1.3180 +            if (ROW_IS_ALIVE (r))
  1.3181 +            {
  1.3182 +                Row [r].shared2.mark = 0 ;
  1.3183 +            }
  1.3184 +        }
  1.3185 +        tag_mark = 1 ;
  1.3186 +    }
  1.3187 +
  1.3188 +    return (tag_mark) ;
  1.3189 +}
  1.3190 +
  1.3191 +
  1.3192 +/* ========================================================================== */
  1.3193 +/* === print_report ========================================================= */
  1.3194 +/* ========================================================================== */
  1.3195 +
  1.3196 +PRIVATE void print_report
  1.3197 +(
  1.3198 +    char *method,
  1.3199 +    Int stats [COLAMD_STATS]
  1.3200 +)
  1.3201 +{
  1.3202 +
  1.3203 +    Int i1, i2, i3 ;
  1.3204 +
  1.3205 +    PRINTF (("\n%s version %d.%d, %s: ", method,
  1.3206 +            COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
  1.3207 +
  1.3208 +    if (!stats)
  1.3209 +    {
  1.3210 +        PRINTF (("No statistics available.\n")) ;
  1.3211 +        return ;
  1.3212 +    }
  1.3213 +
  1.3214 +    i1 = stats [COLAMD_INFO1] ;
  1.3215 +    i2 = stats [COLAMD_INFO2] ;
  1.3216 +    i3 = stats [COLAMD_INFO3] ;
  1.3217 +
  1.3218 +    if (stats [COLAMD_STATUS] >= 0)
  1.3219 +    {
  1.3220 +        PRINTF (("OK.  ")) ;
  1.3221 +    }
  1.3222 +    else
  1.3223 +    {
  1.3224 +        PRINTF (("ERROR.  ")) ;
  1.3225 +    }
  1.3226 +
  1.3227 +    switch (stats [COLAMD_STATUS])
  1.3228 +    {
  1.3229 +
  1.3230 +        case COLAMD_OK_BUT_JUMBLED:
  1.3231 +
  1.3232 +            PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
  1.3233 +
  1.3234 +            PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
  1.3235 +            method, i3)) ;
  1.3236 +
  1.3237 +            PRINTF(("%s: last seen duplicate or out-of-order row index:   %d\n",
  1.3238 +            method, INDEX (i2))) ;
  1.3239 +
  1.3240 +            PRINTF(("%s: last seen in column:                             %d",
  1.3241 +            method, INDEX (i1))) ;
  1.3242 +
  1.3243 +            /* no break - fall through to next case instead */
  1.3244 +
  1.3245 +        case COLAMD_OK:
  1.3246 +
  1.3247 +            PRINTF(("\n")) ;
  1.3248 +
  1.3249 +            PRINTF(("%s: number of dense or empty rows ignored:           %d\n",
  1.3250 +            method, stats [COLAMD_DENSE_ROW])) ;
  1.3251 +
  1.3252 +            PRINTF(("%s: number of dense or empty columns ignored:        %d\n",
  1.3253 +            method, stats [COLAMD_DENSE_COL])) ;
  1.3254 +
  1.3255 +            PRINTF(("%s: number of garbage collections performed:         %d\n",
  1.3256 +            method, stats [COLAMD_DEFRAG_COUNT])) ;
  1.3257 +            break ;
  1.3258 +
  1.3259 +        case COLAMD_ERROR_A_not_present:
  1.3260 +
  1.3261 +            PRINTF(("Array A (row indices of matrix) not present.\n")) ;
  1.3262 +            break ;
  1.3263 +
  1.3264 +        case COLAMD_ERROR_p_not_present:
  1.3265 +
  1.3266 +            PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
  1.3267 +            break ;
  1.3268 +
  1.3269 +        case COLAMD_ERROR_nrow_negative:
  1.3270 +
  1.3271 +            PRINTF(("Invalid number of rows (%d).\n", i1)) ;
  1.3272 +            break ;
  1.3273 +
  1.3274 +        case COLAMD_ERROR_ncol_negative:
  1.3275 +
  1.3276 +            PRINTF(("Invalid number of columns (%d).\n", i1)) ;
  1.3277 +            break ;
  1.3278 +
  1.3279 +        case COLAMD_ERROR_nnz_negative:
  1.3280 +
  1.3281 +            PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
  1.3282 +            break ;
  1.3283 +
  1.3284 +        case COLAMD_ERROR_p0_nonzero:
  1.3285 +
  1.3286 +            PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
  1.3287 +            break ;
  1.3288 +
  1.3289 +        case COLAMD_ERROR_A_too_small:
  1.3290 +
  1.3291 +            PRINTF(("Array A too small.\n")) ;
  1.3292 +            PRINTF(("        Need Alen >= %d, but given only Alen = %d.\n",
  1.3293 +            i1, i2)) ;
  1.3294 +            break ;
  1.3295 +
  1.3296 +        case COLAMD_ERROR_col_length_negative:
  1.3297 +
  1.3298 +            PRINTF
  1.3299 +            (("Column %d has a negative number of nonzero entries (%d).\n",
  1.3300 +            INDEX (i1), i2)) ;
  1.3301 +            break ;
  1.3302 +
  1.3303 +        case COLAMD_ERROR_row_index_out_of_bounds:
  1.3304 +
  1.3305 +            PRINTF
  1.3306 +            (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
  1.3307 +            INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
  1.3308 +            break ;
  1.3309 +
  1.3310 +        case COLAMD_ERROR_out_of_memory:
  1.3311 +
  1.3312 +            PRINTF(("Out of memory.\n")) ;
  1.3313 +            break ;
  1.3314 +
  1.3315 +        /* v2.4: internal-error case deleted */
  1.3316 +    }
  1.3317 +}
  1.3318 +
  1.3319 +
  1.3320 +
  1.3321 +
  1.3322 +/* ========================================================================== */
  1.3323 +/* === colamd debugging routines ============================================ */
  1.3324 +/* ========================================================================== */
  1.3325 +
  1.3326 +/* When debugging is disabled, the remainder of this file is ignored. */
  1.3327 +
  1.3328 +#ifndef NDEBUG
  1.3329 +
  1.3330 +
  1.3331 +/* ========================================================================== */
  1.3332 +/* === debug_structures ===================================================== */
  1.3333 +/* ========================================================================== */
  1.3334 +
  1.3335 +/*
  1.3336 +    At this point, all empty rows and columns are dead.  All live columns
  1.3337 +    are "clean" (containing no dead rows) and simplicial (no supercolumns
  1.3338 +    yet).  Rows may contain dead columns, but all live rows contain at
  1.3339 +    least one live column.
  1.3340 +*/
  1.3341 +
  1.3342 +PRIVATE void debug_structures
  1.3343 +(
  1.3344 +    /* === Parameters ======================================================= */
  1.3345 +
  1.3346 +    Int n_row,
  1.3347 +    Int n_col,
  1.3348 +    Colamd_Row Row [],
  1.3349 +    Colamd_Col Col [],
  1.3350 +    Int A [],
  1.3351 +    Int n_col2
  1.3352 +)
  1.3353 +{
  1.3354 +    /* === Local variables ================================================== */
  1.3355 +
  1.3356 +    Int i ;
  1.3357 +    Int c ;
  1.3358 +    Int *cp ;
  1.3359 +    Int *cp_end ;
  1.3360 +    Int len ;
  1.3361 +    Int score ;
  1.3362 +    Int r ;
  1.3363 +    Int *rp ;
  1.3364 +    Int *rp_end ;
  1.3365 +    Int deg ;
  1.3366 +
  1.3367 +    /* === Check A, Row, and Col ============================================ */
  1.3368 +
  1.3369 +    for (c = 0 ; c < n_col ; c++)
  1.3370 +    {
  1.3371 +        if (COL_IS_ALIVE (c))
  1.3372 +        {
  1.3373 +            len = Col [c].length ;
  1.3374 +            score = Col [c].shared2.score ;
  1.3375 +            DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
  1.3376 +            ASSERT (len > 0) ;
  1.3377 +            ASSERT (score >= 0) ;
  1.3378 +            ASSERT (Col [c].shared1.thickness == 1) ;
  1.3379 +            cp = &A [Col [c].start] ;
  1.3380 +            cp_end = cp + len ;
  1.3381 +            while (cp < cp_end)
  1.3382 +            {
  1.3383 +                r = *cp++ ;
  1.3384 +                ASSERT (ROW_IS_ALIVE (r)) ;
  1.3385 +            }
  1.3386 +        }
  1.3387 +        else
  1.3388 +        {
  1.3389 +            i = Col [c].shared2.order ;
  1.3390 +            ASSERT (i >= n_col2 && i < n_col) ;
  1.3391 +        }
  1.3392 +    }
  1.3393 +
  1.3394 +    for (r = 0 ; r < n_row ; r++)
  1.3395 +    {
  1.3396 +        if (ROW_IS_ALIVE (r))
  1.3397 +        {
  1.3398 +            i = 0 ;
  1.3399 +            len = Row [r].length ;
  1.3400 +            deg = Row [r].shared1.degree ;
  1.3401 +            ASSERT (len > 0) ;
  1.3402 +            ASSERT (deg > 0) ;
  1.3403 +            rp = &A [Row [r].start] ;
  1.3404 +            rp_end = rp + len ;
  1.3405 +            while (rp < rp_end)
  1.3406 +            {
  1.3407 +                c = *rp++ ;
  1.3408 +                if (COL_IS_ALIVE (c))
  1.3409 +                {
  1.3410 +                    i++ ;
  1.3411 +                }
  1.3412 +            }
  1.3413 +            ASSERT (i > 0) ;
  1.3414 +        }
  1.3415 +    }
  1.3416 +}
  1.3417 +
  1.3418 +
  1.3419 +/* ========================================================================== */
  1.3420 +/* === debug_deg_lists ====================================================== */
  1.3421 +/* ========================================================================== */
  1.3422 +
  1.3423 +/*
  1.3424 +    Prints the contents of the degree lists.  Counts the number of columns
  1.3425 +    in the degree list and compares it to the total it should have.  Also
  1.3426 +    checks the row degrees.
  1.3427 +*/
  1.3428 +
  1.3429 +PRIVATE void debug_deg_lists
  1.3430 +(
  1.3431 +    /* === Parameters ======================================================= */
  1.3432 +
  1.3433 +    Int n_row,
  1.3434 +    Int n_col,
  1.3435 +    Colamd_Row Row [],
  1.3436 +    Colamd_Col Col [],
  1.3437 +    Int head [],
  1.3438 +    Int min_score,
  1.3439 +    Int should,
  1.3440 +    Int max_deg
  1.3441 +)
  1.3442 +{
  1.3443 +    /* === Local variables ================================================== */
  1.3444 +
  1.3445 +    Int deg ;
  1.3446 +    Int col ;
  1.3447 +    Int have ;
  1.3448 +    Int row ;
  1.3449 +
  1.3450 +    /* === Check the degree lists =========================================== */
  1.3451 +
  1.3452 +    if (n_col > 10000 && colamd_debug <= 0)
  1.3453 +    {
  1.3454 +        return ;
  1.3455 +    }
  1.3456 +    have = 0 ;
  1.3457 +    DEBUG4 (("Degree lists: %d\n", min_score)) ;
  1.3458 +    for (deg = 0 ; deg <= n_col ; deg++)
  1.3459 +    {
  1.3460 +        col = head [deg] ;
  1.3461 +        if (col == EMPTY)
  1.3462 +        {
  1.3463 +            continue ;
  1.3464 +        }
  1.3465 +        DEBUG4 (("%d:", deg)) ;
  1.3466 +        while (col != EMPTY)
  1.3467 +        {
  1.3468 +            DEBUG4 ((" %d", col)) ;
  1.3469 +            have += Col [col].shared1.thickness ;
  1.3470 +            ASSERT (COL_IS_ALIVE (col)) ;
  1.3471 +            col = Col [col].shared4.degree_next ;
  1.3472 +        }
  1.3473 +        DEBUG4 (("\n")) ;
  1.3474 +    }
  1.3475 +    DEBUG4 (("should %d have %d\n", should, have)) ;
  1.3476 +    ASSERT (should == have) ;
  1.3477 +
  1.3478 +    /* === Check the row degrees ============================================ */
  1.3479 +
  1.3480 +    if (n_row > 10000 && colamd_debug <= 0)
  1.3481 +    {
  1.3482 +        return ;
  1.3483 +    }
  1.3484 +    for (row = 0 ; row < n_row ; row++)
  1.3485 +    {
  1.3486 +        if (ROW_IS_ALIVE (row))
  1.3487 +        {
  1.3488 +            ASSERT (Row [row].shared1.degree <= max_deg) ;
  1.3489 +        }
  1.3490 +    }
  1.3491 +}
  1.3492 +
  1.3493 +
  1.3494 +/* ========================================================================== */
  1.3495 +/* === debug_mark =========================================================== */
  1.3496 +/* ========================================================================== */
  1.3497 +
  1.3498 +/*
  1.3499 +    Ensures that the tag_mark is less that the maximum and also ensures that
  1.3500 +    each entry in the mark array is less than the tag mark.
  1.3501 +*/
  1.3502 +
  1.3503 +PRIVATE void debug_mark
  1.3504 +(
  1.3505 +    /* === Parameters ======================================================= */
  1.3506 +
  1.3507 +    Int n_row,
  1.3508 +    Colamd_Row Row [],
  1.3509 +    Int tag_mark,
  1.3510 +    Int max_mark
  1.3511 +)
  1.3512 +{
  1.3513 +    /* === Local variables ================================================== */
  1.3514 +
  1.3515 +    Int r ;
  1.3516 +
  1.3517 +    /* === Check the Row marks ============================================== */
  1.3518 +
  1.3519 +    ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
  1.3520 +    if (n_row > 10000 && colamd_debug <= 0)
  1.3521 +    {
  1.3522 +        return ;
  1.3523 +    }
  1.3524 +    for (r = 0 ; r < n_row ; r++)
  1.3525 +    {
  1.3526 +        ASSERT (Row [r].shared2.mark < tag_mark) ;
  1.3527 +    }
  1.3528 +}
  1.3529 +
  1.3530 +
  1.3531 +/* ========================================================================== */
  1.3532 +/* === debug_matrix ========================================================= */
  1.3533 +/* ========================================================================== */
  1.3534 +
  1.3535 +/*
  1.3536 +    Prints out the contents of the columns and the rows.
  1.3537 +*/
  1.3538 +
  1.3539 +PRIVATE void debug_matrix
  1.3540 +(
  1.3541 +    /* === Parameters ======================================================= */
  1.3542 +
  1.3543 +    Int n_row,
  1.3544 +    Int n_col,
  1.3545 +    Colamd_Row Row [],
  1.3546 +    Colamd_Col Col [],
  1.3547 +    Int A []
  1.3548 +)
  1.3549 +{
  1.3550 +    /* === Local variables ================================================== */
  1.3551 +
  1.3552 +    Int r ;
  1.3553 +    Int c ;
  1.3554 +    Int *rp ;
  1.3555 +    Int *rp_end ;
  1.3556 +    Int *cp ;
  1.3557 +    Int *cp_end ;
  1.3558 +
  1.3559 +    /* === Dump the rows and columns of the matrix ========================== */
  1.3560 +
  1.3561 +    if (colamd_debug < 3)
  1.3562 +    {
  1.3563 +        return ;
  1.3564 +    }
  1.3565 +    DEBUG3 (("DUMP MATRIX:\n")) ;
  1.3566 +    for (r = 0 ; r < n_row ; r++)
  1.3567 +    {
  1.3568 +        DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
  1.3569 +        if (ROW_IS_DEAD (r))
  1.3570 +        {
  1.3571 +            continue ;
  1.3572 +        }
  1.3573 +        DEBUG3 (("start %d length %d degree %d\n",
  1.3574 +                Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
  1.3575 +        rp = &A [Row [r].start] ;
  1.3576 +        rp_end = rp + Row [r].length ;
  1.3577 +        while (rp < rp_end)
  1.3578 +        {
  1.3579 +            c = *rp++ ;
  1.3580 +            DEBUG4 (("  %d col %d\n", COL_IS_ALIVE (c), c)) ;
  1.3581 +        }
  1.3582 +    }
  1.3583 +
  1.3584 +    for (c = 0 ; c < n_col ; c++)
  1.3585 +    {
  1.3586 +        DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
  1.3587 +        if (COL_IS_DEAD (c))
  1.3588 +        {
  1.3589 +            continue ;
  1.3590 +        }
  1.3591 +        DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
  1.3592 +                Col [c].start, Col [c].length,
  1.3593 +                Col [c].shared1.thickness, Col [c].shared2.score)) ;
  1.3594 +        cp = &A [Col [c].start] ;
  1.3595 +        cp_end = cp + Col [c].length ;
  1.3596 +        while (cp < cp_end)
  1.3597 +        {
  1.3598 +            r = *cp++ ;
  1.3599 +            DEBUG4 (("  %d row %d\n", ROW_IS_ALIVE (r), r)) ;
  1.3600 +        }
  1.3601 +    }
  1.3602 +}
  1.3603 +
  1.3604 +PRIVATE void colamd_get_debug
  1.3605 +(
  1.3606 +    char *method
  1.3607 +)
  1.3608 +{
  1.3609 +    FILE *f ;
  1.3610 +    colamd_debug = 0 ;          /* no debug printing */
  1.3611 +    f = fopen ("debug", "r") ;
  1.3612 +    if (f == (FILE *) NULL)
  1.3613 +    {
  1.3614 +        colamd_debug = 0 ;
  1.3615 +    }
  1.3616 +    else
  1.3617 +    {
  1.3618 +        fscanf (f, "%d", &colamd_debug) ;
  1.3619 +        fclose (f) ;
  1.3620 +    }
  1.3621 +    DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
  1.3622 +        method, colamd_debug)) ;
  1.3623 +}
  1.3624 +
  1.3625 +#endif /* NDEBUG */