Amesos Package Browser (Single Doxygen Collection) Development
amesos_colamd_l.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
00003 /* ========================================================================== */
00004 
00005 /* COLAMD / SYMAMD
00006 
00007     colamd:  an approximate minimum degree column ordering algorithm,
00008       for LU factorization of symmetric or unsymmetric matrices,
00009   QR factorization, least squares, interior point methods for
00010   linear programming problems, and other related problems.
00011 
00012     symamd:  an approximate minimum degree ordering algorithm for Cholesky
00013       factorization of symmetric matrices.
00014 
00015     Purpose:
00016 
00017   Colamd computes a permutation Q such that the Cholesky factorization of
00018   (AQ)'(AQ) has less fill-in and requires fewer floating point operations
00019   than A'A.  This also provides a good ordering for sparse partial
00020   pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
00021   factorization, and P is computed during numerical factorization via
00022   conventional partial pivoting with row interchanges.  Colamd is the
00023   column ordering method used in SuperLU, part of the ScaLAPACK library.
00024   It is also available as built-in function in MATLAB Version 6,
00025   available from MathWorks, Inc. (http://www.mathworks.com).  This
00026   routine can be used in place of colmmd in MATLAB.
00027 
00028       Symamd computes a permutation P of a symmetric matrix A such that the
00029   Cholesky factorization of PAP' has less fill-in and requires fewer
00030   floating point operations than A.  Symamd constructs a matrix M such
00031   that M'M has the same nonzero pattern of A, and then orders the columns
00032   of M using colmmd.  The column ordering of M is then returned as the
00033   row and column ordering P of A. 
00034 
00035     Authors:
00036 
00037   The authors of the code itself are Stefan I. Larimore and Timothy A.
00038   Davis (davis at cise.ufl.edu), University of Florida.  The algorithm was
00039   developed in collaboration with John Gilbert, Xerox PARC, and Esmond
00040   Ng, Oak Ridge National Laboratory.
00041 
00042     Acknowledgements:
00043 
00044   This work was supported by the National Science Foundation, under
00045   grants DMS-9504974 and DMS-9803599.
00046 
00047     Copyright and License:
00048 
00049   Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
00050   COLAMD is also available under alternate licenses, contact T. Davis
00051   for details.
00052 
00053   This library is free software; you can redistribute it and/or
00054   modify it under the terms of the GNU Lesser General Public
00055   License as published by the Free Software Foundation; either
00056   version 2.1 of the License, or (at your option) any later version.
00057 
00058   This library is distributed in the hope that it will be useful,
00059   but WITHOUT ANY WARRANTY; without even the implied warranty of
00060   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00061   Lesser General Public License for more details.
00062 
00063   You should have received a copy of the GNU Lesser General Public
00064   License along with this library; if not, write to the Free Software
00065   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
00066   USA
00067 
00068   Permission is hereby granted to use or copy this program under the
00069   terms of the GNU LGPL, provided that the Copyright, this License,
00070   and the Availability of the original version is retained on all copies.
00071   User documentation of any code that uses this code or any modified
00072   version of this code must cite the Copyright, this License, the
00073   Availability note, and "Used by permission." Permission to modify
00074   the code and to distribute modified code is granted, provided the
00075   Copyright, this License, and the Availability note are retained,
00076   and a notice that the code was modified is included.
00077 
00078     Availability:
00079 
00080   The colamd/symamd library is available at
00081 
00082       http://www.cise.ufl.edu/research/sparse/colamd/
00083 
00084   This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
00085   file.  It requires the colamd.h file.  It is required by the colamdmex.c
00086   and symamdmex.c files, for the MATLAB interface to colamd and symamd.
00087   Appears as ACM Algorithm 836.
00088 
00089     See the ChangeLog file for changes since Version 1.0.
00090 
00091     References:
00092 
00093   T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
00094   minimum degree ordering algorithm, ACM Transactions on Mathematical
00095   Software, vol. 30, no. 3., pp. 353-376, 2004.
00096 
00097   T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
00098   an approximate column minimum degree ordering algorithm, ACM
00099   Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
00100   2004.
00101 
00102 */
00103 
00104 /* ========================================================================== */
00105 /* === Description of user-callable routines ================================ */
00106 /* ========================================================================== */
00107 
00108 /* COLAMD includes both int and UF_long versions of all its routines.  The
00109  * description below is for the int version.  For UF_long, all int arguments
00110  * become UF_long.  UF_long is normally defined as long, except for WIN64.
00111 
00112     ----------------------------------------------------------------------------
00113     colamd_recommended:
00114     ----------------------------------------------------------------------------
00115 
00116   C syntax:
00117 
00118       #include "colamd.h"
00119       size_t colamd_recommended (int nnz, int n_row, int n_col) ;
00120       size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
00121     UF_long n_col) ;
00122 
00123   Purpose:
00124 
00125       Returns recommended value of Alen for use by colamd.  Returns 0
00126       if any input argument is negative.  The use of this routine
00127       is optional.  Not needed for symamd, which dynamically allocates
00128       its own memory.
00129 
00130       Note that in v2.4 and earlier, these routines returned int or long.
00131       They now return a value of type size_t.
00132 
00133   Arguments (all input arguments):
00134 
00135       int nnz ;   Number of nonzeros in the matrix A.  This must
00136         be the same value as p [n_col] in the call to
00137         colamd - otherwise you will get a wrong value
00138         of the recommended memory to use.
00139 
00140       int n_row ;   Number of rows in the matrix A.
00141 
00142       int n_col ;   Number of columns in the matrix A.
00143 
00144     ----------------------------------------------------------------------------
00145     colamd_set_defaults:
00146     ----------------------------------------------------------------------------
00147 
00148   C syntax:
00149 
00150       #include "colamd.h"
00151       colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
00152       colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
00153 
00154   Purpose:
00155 
00156       Sets the default parameters.  The use of this routine is optional.
00157 
00158   Arguments:
00159 
00160       double knobs [COLAMD_KNOBS] ; Output only.
00161 
00162     NOTE: the meaning of the dense row/col knobs has changed in v2.4
00163 
00164     knobs [0] and knobs [1] control dense row and col detection:
00165 
00166     Colamd: rows with more than
00167     max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
00168     entries are removed prior to ordering.  Columns with more than
00169     max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
00170     entries are removed prior to
00171     ordering, and placed last in the output column ordering. 
00172 
00173     Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
00174     Rows and columns with more than
00175     max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
00176     entries are removed prior to ordering, and placed last in the
00177     output ordering.
00178 
00179     COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
00180     respectively, in colamd.h.  Default values of these two knobs
00181     are both 10.  Currently, only knobs [0] and knobs [1] are
00182     used, but future versions may use more knobs.  If so, they will
00183     be properly set to their defaults by the future version of
00184     colamd_set_defaults, so that the code that calls colamd will
00185     not need to change, assuming that you either use
00186     colamd_set_defaults, or pass a (double *) NULL pointer as the
00187     knobs array to colamd or symamd.
00188 
00189       knobs [2]: aggressive absorption
00190 
00191           knobs [COLAMD_AGGRESSIVE] controls whether or not to do
00192           aggressive absorption during the ordering.  Default is TRUE.
00193 
00194 
00195     ----------------------------------------------------------------------------
00196     colamd:
00197     ----------------------------------------------------------------------------
00198 
00199   C syntax:
00200 
00201       #include "colamd.h"
00202       int colamd (int n_row, int n_col, int Alen, int *A, int *p,
00203         double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
00204       UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
00205     UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
00206     UF_long stats [COLAMD_STATS]) ;
00207 
00208   Purpose:
00209 
00210       Computes a column ordering (Q) of A such that P(AQ)=LU or
00211       (AQ)'AQ=LL' have less fill-in and require fewer floating point
00212       operations than factorizing the unpermuted matrix A or A'A,
00213       respectively.
00214       
00215   Returns:
00216 
00217       TRUE (1) if successful, FALSE (0) otherwise.
00218 
00219   Arguments:
00220 
00221       int n_row ;   Input argument.
00222 
00223     Number of rows in the matrix A.
00224     Restriction:  n_row >= 0.
00225     Colamd returns FALSE if n_row is negative.
00226 
00227       int n_col ;   Input argument.
00228 
00229     Number of columns in the matrix A.
00230     Restriction:  n_col >= 0.
00231     Colamd returns FALSE if n_col is negative.
00232 
00233       int Alen ;    Input argument.
00234 
00235     Restriction (see note):
00236     Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
00237     Colamd returns FALSE if these conditions are not met.
00238 
00239     Note:  this restriction makes an modest assumption regarding
00240     the size of the two typedef's structures in colamd.h.
00241     We do, however, guarantee that
00242 
00243       Alen >= colamd_recommended (nnz, n_row, n_col)
00244 
00245     will be sufficient.  Note: the macro version does not check
00246     for integer overflow, and thus is not recommended.  Use
00247     the colamd_recommended routine instead.
00248 
00249       int A [Alen] ;  Input argument, undefined on output.
00250 
00251     A is an integer array of size Alen.  Alen must be at least as
00252     large as the bare minimum value given above, but this is very
00253     low, and can result in excessive run time.  For best
00254     performance, we recommend that Alen be greater than or equal to
00255     colamd_recommended (nnz, n_row, n_col), which adds
00256     nnz/5 to the bare minimum value given above.
00257 
00258     On input, the row indices of the entries in column c of the
00259     matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
00260     in a given column c need not be in ascending order, and
00261     duplicate row indices may be be present.  However, colamd will
00262     work a little faster if both of these conditions are met
00263     (Colamd puts the matrix into this format, if it finds that the
00264     the conditions are not met).
00265 
00266     The matrix is 0-based.  That is, rows are in the range 0 to
00267     n_row-1, and columns are in the range 0 to n_col-1.  Colamd
00268     returns FALSE if any row index is out of range.
00269 
00270     The contents of A are modified during ordering, and are
00271     undefined on output.
00272 
00273       int p [n_col+1] ; Both input and output argument.
00274 
00275     p is an integer array of size n_col+1.  On input, it holds the
00276     "pointers" for the column form of the matrix A.  Column c of
00277     the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
00278     entry, p [0], must be zero, and p [c] <= p [c+1] must hold
00279     for all c in the range 0 to n_col-1.  The value p [n_col] is
00280     thus the total number of entries in the pattern of the matrix A.
00281     Colamd returns FALSE if these conditions are not met.
00282 
00283     On output, if colamd returns TRUE, the array p holds the column
00284     permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
00285     the first column index in the new ordering, and p [n_col-1] is
00286     the last.  That is, p [k] = j means that column j of A is the
00287     kth pivot column, in AQ, where k is in the range 0 to n_col-1
00288     (p [0] = j means that column j of A is the first column in AQ).
00289 
00290     If colamd returns FALSE, then no permutation is returned, and
00291     p is undefined on output.
00292 
00293       double knobs [COLAMD_KNOBS] ; Input argument.
00294 
00295     See colamd_set_defaults for a description.
00296 
00297       int stats [COLAMD_STATS] ;    Output argument.
00298 
00299     Statistics on the ordering, and error status.
00300     See colamd.h for related definitions.
00301     Colamd returns FALSE if stats is not present.
00302 
00303     stats [0]:  number of dense or empty rows ignored.
00304 
00305     stats [1]:  number of dense or empty columns ignored (and
00306         ordered last in the output permutation p)
00307         Note that a row can become "empty" if it
00308         contains only "dense" and/or "empty" columns,
00309         and similarly a column can become "empty" if it
00310         only contains "dense" and/or "empty" rows.
00311 
00312     stats [2]:  number of garbage collections performed.
00313         This can be excessively high if Alen is close
00314         to the minimum required value.
00315 
00316     stats [3]:  status code.  < 0 is an error code.
00317           > 1 is a warning or notice.
00318 
00319       0 OK.  Each column of the input matrix contained
00320         row indices in increasing order, with no
00321         duplicates.
00322 
00323       1 OK, but columns of input matrix were jumbled
00324         (unsorted columns or duplicate entries).  Colamd
00325         had to do some extra work to sort the matrix
00326         first and remove duplicate entries, but it
00327         still was able to return a valid permutation
00328         (return value of colamd was TRUE).
00329 
00330           stats [4]: highest numbered column that
00331             is unsorted or has duplicate
00332             entries.
00333           stats [5]: last seen duplicate or
00334             unsorted row index.
00335           stats [6]: number of duplicate or
00336             unsorted row indices.
00337 
00338       -1  A is a null pointer
00339 
00340       -2  p is a null pointer
00341 
00342       -3  n_row is negative
00343 
00344           stats [4]: n_row
00345 
00346       -4  n_col is negative
00347 
00348           stats [4]: n_col
00349 
00350       -5  number of nonzeros in matrix is negative
00351 
00352           stats [4]: number of nonzeros, p [n_col]
00353 
00354       -6  p [0] is nonzero
00355 
00356           stats [4]: p [0]
00357 
00358       -7  A is too small
00359 
00360           stats [4]: required size
00361           stats [5]: actual size (Alen)
00362 
00363       -8  a column has a negative number of entries
00364 
00365           stats [4]: column with < 0 entries
00366           stats [5]: number of entries in col
00367 
00368       -9  a row index is out of bounds
00369 
00370           stats [4]: column with bad row index
00371           stats [5]: bad row index
00372           stats [6]: n_row, # of rows of matrx
00373 
00374       -10 (unused; see symamd.c)
00375 
00376       -999  (unused; see symamd.c)
00377 
00378     Future versions may return more statistics in the stats array.
00379 
00380   Example:
00381   
00382       See http://www.cise.ufl.edu/research/sparse/colamd/example.c
00383       for a complete example.
00384 
00385       To order the columns of a 5-by-4 matrix with 11 nonzero entries in
00386       the following nonzero pattern
00387 
00388         x 0 x 0
00389     x 0 x x
00390     0 x x 0
00391     0 0 x x
00392     x x 0 0
00393 
00394       with default knobs and no output statistics, do the following:
00395 
00396     #include "colamd.h"
00397     #define ALEN 100
00398     int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
00399     int p [ ] = {0, 3, 5, 9, 11} ;
00400     int stats [COLAMD_STATS] ;
00401     colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
00402 
00403       The permutation is returned in the array p, and A is destroyed.
00404 
00405     ----------------------------------------------------------------------------
00406     symamd:
00407     ----------------------------------------------------------------------------
00408 
00409   C syntax:
00410 
00411       #include "colamd.h"
00412       int symamd (int n, int *A, int *p, int *perm,
00413         double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
00414     void (*allocate) (size_t, size_t), void (*release) (void *)) ;
00415       UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
00416         double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
00417     void (*allocate) (size_t, size_t), void (*release) (void *)) ;
00418 
00419   Purpose:
00420 
00421           The symamd routine computes an ordering P of a symmetric sparse
00422       matrix A such that the Cholesky factorization PAP' = LL' remains
00423       sparse.  It is based on a column ordering of a matrix M constructed
00424       so that the nonzero pattern of M'M is the same as A.  The matrix A
00425       is assumed to be symmetric; only the strictly lower triangular part
00426       is accessed.  You must pass your selected memory allocator (usually
00427       calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
00428       memory for the temporary matrix M.
00429 
00430   Returns:
00431 
00432       TRUE (1) if successful, FALSE (0) otherwise.
00433 
00434   Arguments:
00435 
00436       int n ;   Input argument.
00437 
00438         Number of rows and columns in the symmetrix matrix A.
00439     Restriction:  n >= 0.
00440     Symamd returns FALSE if n is negative.
00441 
00442       int A [nnz] ; Input argument.
00443 
00444         A is an integer array of size nnz, where nnz = p [n].
00445     
00446     The row indices of the entries in column c of the matrix are
00447     held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
00448     given column c need not be in ascending order, and duplicate
00449     row indices may be present.  However, symamd will run faster
00450     if the columns are in sorted order with no duplicate entries. 
00451 
00452     The matrix is 0-based.  That is, rows are in the range 0 to
00453     n-1, and columns are in the range 0 to n-1.  Symamd
00454     returns FALSE if any row index is out of range.
00455 
00456     The contents of A are not modified.
00457 
00458       int p [n+1] ;     Input argument.
00459 
00460     p is an integer array of size n+1.  On input, it holds the
00461     "pointers" for the column form of the matrix A.  Column c of
00462     the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
00463     entry, p [0], must be zero, and p [c] <= p [c+1] must hold
00464     for all c in the range 0 to n-1.  The value p [n] is
00465     thus the total number of entries in the pattern of the matrix A.
00466     Symamd returns FALSE if these conditions are not met.
00467 
00468     The contents of p are not modified.
00469 
00470       int perm [n+1] ;    Output argument.
00471 
00472     On output, if symamd returns TRUE, the array perm holds the
00473     permutation P, where perm [0] is the first index in the new
00474     ordering, and perm [n-1] is the last.  That is, perm [k] = j
00475     means that row and column j of A is the kth column in PAP',
00476     where k is in the range 0 to n-1 (perm [0] = j means
00477     that row and column j of A are the first row and column in
00478     PAP').  The array is used as a workspace during the ordering,
00479     which is why it must be of length n+1, not just n.
00480 
00481       double knobs [COLAMD_KNOBS] ; Input argument.
00482 
00483     See colamd_set_defaults for a description.
00484 
00485       int stats [COLAMD_STATS] ;    Output argument.
00486 
00487     Statistics on the ordering, and error status.
00488     See colamd.h for related definitions.
00489     Symamd returns FALSE if stats is not present.
00490 
00491     stats [0]:  number of dense or empty row and columns ignored
00492         (and ordered last in the output permutation 
00493         perm).  Note that a row/column can become
00494         "empty" if it contains only "dense" and/or
00495         "empty" columns/rows.
00496 
00497     stats [1]:  (same as stats [0])
00498 
00499     stats [2]:  number of garbage collections performed.
00500 
00501     stats [3]:  status code.  < 0 is an error code.
00502           > 1 is a warning or notice.
00503 
00504       0 OK.  Each column of the input matrix contained
00505         row indices in increasing order, with no
00506         duplicates.
00507 
00508       1 OK, but columns of input matrix were jumbled
00509         (unsorted columns or duplicate entries).  Symamd
00510         had to do some extra work to sort the matrix
00511         first and remove duplicate entries, but it
00512         still was able to return a valid permutation
00513         (return value of symamd was TRUE).
00514 
00515           stats [4]: highest numbered column that
00516             is unsorted or has duplicate
00517             entries.
00518           stats [5]: last seen duplicate or
00519             unsorted row index.
00520           stats [6]: number of duplicate or
00521             unsorted row indices.
00522 
00523       -1  A is a null pointer
00524 
00525       -2  p is a null pointer
00526 
00527       -3  (unused, see colamd.c)
00528 
00529       -4  n is negative
00530 
00531           stats [4]: n
00532 
00533       -5  number of nonzeros in matrix is negative
00534 
00535           stats [4]: # of nonzeros (p [n]).
00536 
00537       -6  p [0] is nonzero
00538 
00539           stats [4]: p [0]
00540 
00541       -7  (unused)
00542 
00543       -8  a column has a negative number of entries
00544 
00545           stats [4]: column with < 0 entries
00546           stats [5]: number of entries in col
00547 
00548       -9  a row index is out of bounds
00549 
00550           stats [4]: column with bad row index
00551           stats [5]: bad row index
00552           stats [6]: n_row, # of rows of matrx
00553 
00554       -10 out of memory (unable to allocate temporary
00555         workspace for M or count arrays using the
00556         "allocate" routine passed into symamd).
00557 
00558     Future versions may return more statistics in the stats array.
00559 
00560       void * (*allocate) (size_t, size_t)
00561 
00562         A pointer to a function providing memory allocation.  The
00563     allocated memory must be returned initialized to zero.  For a
00564     C application, this argument should normally be a pointer to
00565     calloc.  For a MATLAB mexFunction, the routine mxCalloc is
00566     passed instead.
00567 
00568       void (*release) (size_t, size_t)
00569 
00570         A pointer to a function that frees memory allocated by the
00571     memory allocation routine above.  For a C application, this
00572     argument should normally be a pointer to free.  For a MATLAB
00573     mexFunction, the routine mxFree is passed instead.
00574 
00575 
00576     ----------------------------------------------------------------------------
00577     colamd_report:
00578     ----------------------------------------------------------------------------
00579 
00580   C syntax:
00581 
00582       #include "colamd.h"
00583       colamd_report (int stats [COLAMD_STATS]) ;
00584       colamd_l_report (UF_long stats [COLAMD_STATS]) ;
00585 
00586   Purpose:
00587 
00588       Prints the error status and statistics recorded in the stats
00589       array on the standard error output (for a standard C routine)
00590       or on the MATLAB output (for a mexFunction).
00591 
00592   Arguments:
00593 
00594       int stats [COLAMD_STATS] ;  Input only.  Statistics from colamd.
00595 
00596 
00597     ----------------------------------------------------------------------------
00598     symamd_report:
00599     ----------------------------------------------------------------------------
00600 
00601   C syntax:
00602 
00603       #include "colamd.h"
00604       symamd_report (int stats [COLAMD_STATS]) ;
00605       symamd_l_report (UF_long stats [COLAMD_STATS]) ;
00606 
00607   Purpose:
00608 
00609       Prints the error status and statistics recorded in the stats
00610       array on the standard error output (for a standard C routine)
00611       or on the MATLAB output (for a mexFunction).
00612 
00613   Arguments:
00614 
00615       int stats [COLAMD_STATS] ;  Input only.  Statistics from symamd.
00616 
00617 
00618 */
00619 
00620 /* ========================================================================== */
00621 /* === Scaffolding code definitions  ======================================== */
00622 /* ========================================================================== */
00623 
00624 /* Ensure that debugging is turned off: */
00625 #ifndef NDEBUG
00626 #define NDEBUG
00627 #endif
00628 
00629 /* turn on debugging by uncommenting the following line
00630  #undef NDEBUG
00631 */
00632 
00633 /*
00634    Our "scaffolding code" philosophy:  In our opinion, well-written library
00635    code should keep its "debugging" code, and just normally have it turned off
00636    by the compiler so as not to interfere with performance.  This serves
00637    several purposes:
00638 
00639    (1) assertions act as comments to the reader, telling you what the code
00640   expects at that point.  All assertions will always be true (unless
00641   there really is a bug, of course).
00642 
00643    (2) leaving in the scaffolding code assists anyone who would like to modify
00644   the code, or understand the algorithm (by reading the debugging output,
00645   one can get a glimpse into what the code is doing).
00646 
00647    (3) (gasp!) for actually finding bugs.  This code has been heavily tested
00648   and "should" be fully functional and bug-free ... but you never know...
00649 
00650     The code will become outrageously slow when debugging is
00651     enabled.  To control the level of debugging output, set an environment
00652     variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
00653     you should see the following message on the standard output:
00654 
00655       colamd: debug version, D = 1 (THIS WILL BE SLOW!)
00656 
00657     or a similar message for symamd.  If you don't, then debugging has not
00658     been enabled.
00659 
00660 */
00661 
00662 /* ========================================================================== */
00663 /* === Include files ======================================================== */
00664 /* ========================================================================== */
00665 
00666 /* This file should make the long int version of COLAMD */
00667 #define DLONG 1
00668 
00669 #include "amesos_colamd.h"
00670 #include <limits.h>
00671 #include <math.h>
00672 
00673 #ifdef MATLAB_MEX_FILE
00674 #include "mex.h"
00675 #include "matrix.h"
00676 #endif /* MATLAB_MEX_FILE */
00677 
00678 #if !defined (NPRINT) || !defined (NDEBUG)
00679 #include <stdio.h>
00680 #endif
00681 
00682 #ifndef NULL
00683 #define NULL ((void *) 0)
00684 #endif
00685 
00686 /* ========================================================================== */
00687 /* === int or UF_long ======================================================= */
00688 /* ========================================================================== */
00689 
00690 /* define UF_long */
00691 #include "amesos_UFconfig.h"
00692 
00693 #ifdef DLONG
00694 
00695 #define Int UF_long
00696 #define ID  UF_long_id
00697 #define Int_MAX UF_long_max
00698 
00699 #define COLAMD_recommended amesos_colamd_l_recommended
00700 #define COLAMD_set_defaults amesos_colamd_l_set_defaults
00701 #define COLAMD_MAIN amesos_colamd_l
00702 #define SYMAMD_MAIN amesos_symamd_l
00703 #define COLAMD_report amesos_colamd_l_report
00704 #define SYMAMD_report amesos_symamd_l_report
00705 
00706 #else
00707 
00708 #define Int int
00709 #define ID "%d"
00710 #define Int_MAX INT_MAX
00711 
00712 #define COLAMD_recommended amesos_colamd_recommended
00713 #define COLAMD_set_defaults amesos_colamd_set_defaults
00714 #define COLAMD_MAIN amesos_colamd
00715 #define SYMAMD_MAIN amesos_symamd
00716 #define COLAMD_report amesos_colamd_report
00717 #define SYMAMD_report amesos_symamd_report
00718 
00719 #endif
00720 
00721 /* ========================================================================== */
00722 /* === Row and Column structures ============================================ */
00723 /* ========================================================================== */
00724 
00725 /* User code that makes use of the colamd/symamd routines need not directly */
00726 /* reference these structures.  They are used only for colamd_recommended. */
00727 
00728 typedef struct Colamd_Col_struct
00729 {
00730     Int start ;   /* index for A of first row in this column, or DEAD */
00731       /* if column is dead */
00732     Int length ;  /* number of rows in this column */
00733     union
00734     {
00735   Int thickness ; /* number of original columns represented by this */
00736       /* col, if the column is alive */
00737   Int parent ;  /* parent in parent tree super-column structure, if */
00738       /* the column is dead */
00739     } shared1 ;
00740     union
00741     {
00742   Int score ; /* the score used to maintain heap, if col is alive */
00743   Int order ; /* pivot ordering of this column, if col is dead */
00744     } shared2 ;
00745     union
00746     {
00747   Int headhash ;  /* head of a hash bucket, if col is at the head of */
00748       /* a degree list */
00749   Int hash ;  /* hash value, if col is not in a degree list */
00750   Int prev ;  /* previous column in degree list, if col is in a */
00751       /* degree list (but not at the head of a degree list) */
00752     } shared3 ;
00753     union
00754     {
00755   Int degree_next ; /* next column, if col is in a degree list */
00756   Int hash_next ;   /* next column, if col is in a hash list */
00757     } shared4 ;
00758 
00759 } Colamd_Col ;
00760 
00761 typedef struct Colamd_Row_struct
00762 {
00763     Int start ;   /* index for A of first col in this row */
00764     Int length ;  /* number of principal columns in this row */
00765     union
00766     {
00767   Int degree ;  /* number of principal & non-principal columns in row */
00768   Int p ;   /* used as a row pointer in init_rows_cols () */
00769     } shared1 ;
00770     union
00771     {
00772   Int mark ;  /* for computing set differences and marking dead rows*/
00773   Int first_column ;/* first column in row (used in garbage collection) */
00774     } shared2 ;
00775 
00776 } Colamd_Row ;
00777 
00778 /* ========================================================================== */
00779 /* === Definitions ========================================================== */
00780 /* ========================================================================== */
00781 
00782 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
00783 #define PUBLIC
00784 #define PRIVATE static
00785 
00786 #define DENSE_DEGREE(alpha,n) \
00787     ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
00788 
00789 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
00790 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
00791 
00792 #define ONES_COMPLEMENT(r) (-(r)-1)
00793 
00794 /* -------------------------------------------------------------------------- */
00795 /* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
00796 /* -------------------------------------------------------------------------- */
00797 
00798 #ifndef TRUE
00799 #define TRUE (1)
00800 #endif
00801 
00802 #ifndef FALSE
00803 #define FALSE (0)
00804 #endif
00805 
00806 /* -------------------------------------------------------------------------- */
00807 
00808 #define EMPTY (-1)
00809 
00810 /* Row and column status */
00811 #define ALIVE (0)
00812 #define DEAD  (-1)
00813 
00814 /* Column status */
00815 #define DEAD_PRINCIPAL    (-1)
00816 #define DEAD_NON_PRINCIPAL  (-2)
00817 
00818 /* Macros for row and column status update and checking. */
00819 #define ROW_IS_DEAD(r)      ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
00820 #define ROW_IS_MARKED_DEAD(row_mark)  (row_mark < ALIVE)
00821 #define ROW_IS_ALIVE(r)     (Row [r].shared2.mark >= ALIVE)
00822 #define COL_IS_DEAD(c)      (Col [c].start < ALIVE)
00823 #define COL_IS_ALIVE(c)     (Col [c].start >= ALIVE)
00824 #define COL_IS_DEAD_PRINCIPAL(c)  (Col [c].start == DEAD_PRINCIPAL)
00825 #define KILL_ROW(r)     { Row [r].shared2.mark = DEAD ; }
00826 #define KILL_PRINCIPAL_COL(c)   { Col [c].start = DEAD_PRINCIPAL ; }
00827 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
00828 
00829 /* ========================================================================== */
00830 /* === Colamd reporting mechanism =========================================== */
00831 /* ========================================================================== */
00832 
00833 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
00834 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
00835 #define INDEX(i) ((i)+1)
00836 #else
00837 /* In C, matrices are 0-based and indices are reported as such in *_report */
00838 #define INDEX(i) (i)
00839 #endif
00840 
00841 /* All output goes through the PRINTF macro.  */
00842 #define PRINTF(params) { if (amesos_colamd_printf != NULL) (void) amesos_colamd_printf params ; }
00843 
00844 /* ========================================================================== */
00845 /* === Prototypes of PRIVATE routines ======================================= */
00846 /* ========================================================================== */
00847 
00848 PRIVATE Int init_rows_cols
00849 (
00850     Int n_row,
00851     Int n_col,
00852     Colamd_Row Row [],
00853     Colamd_Col Col [],
00854     Int A [],
00855     Int p [],
00856     Int stats [COLAMD_STATS]
00857 ) ;
00858 
00859 PRIVATE void init_scoring
00860 (
00861     Int n_row,
00862     Int n_col,
00863     Colamd_Row Row [],
00864     Colamd_Col Col [],
00865     Int A [],
00866     Int head [],
00867     double knobs [COLAMD_KNOBS],
00868     Int *p_n_row2,
00869     Int *p_n_col2,
00870     Int *p_max_deg
00871 ) ;
00872 
00873 PRIVATE Int find_ordering
00874 (
00875     Int n_row,
00876     Int n_col,
00877     Int Alen,
00878     Colamd_Row Row [],
00879     Colamd_Col Col [],
00880     Int A [],
00881     Int head [],
00882     Int n_col2,
00883     Int max_deg,
00884     Int pfree,
00885     Int aggressive
00886 ) ;
00887 
00888 PRIVATE void order_children
00889 (
00890     Int n_col,
00891     Colamd_Col Col [],
00892     Int p []
00893 ) ;
00894 
00895 PRIVATE void detect_super_cols
00896 (
00897 
00898 #ifndef NDEBUG
00899     Int n_col,
00900     Colamd_Row Row [],
00901 #endif /* NDEBUG */
00902 
00903     Colamd_Col Col [],
00904     Int A [],
00905     Int head [],
00906     Int row_start,
00907     Int row_length
00908 ) ;
00909 
00910 PRIVATE Int garbage_collection
00911 (
00912     Int n_row,
00913     Int n_col,
00914     Colamd_Row Row [],
00915     Colamd_Col Col [],
00916     Int A [],
00917     Int *pfree
00918 ) ;
00919 
00920 PRIVATE Int clear_mark
00921 (
00922     Int tag_mark,
00923     Int max_mark,
00924     Int n_row,
00925     Colamd_Row Row []
00926 ) ;
00927 
00928 PRIVATE void print_report
00929 (
00930     char *method,
00931     Int stats [COLAMD_STATS]
00932 ) ;
00933 
00934 /* ========================================================================== */
00935 /* === Debugging prototypes and definitions ================================= */
00936 /* ========================================================================== */
00937 
00938 #ifndef NDEBUG
00939 
00940 #include <assert.h>
00941 
00942 /* colamd_debug is the *ONLY* global variable, and is only */
00943 /* present when debugging */
00944 
00945 PRIVATE Int colamd_debug = 0 ;  /* debug print level */
00946 
00947 #define DEBUG0(params) { PRINTF (params) ; }
00948 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
00949 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
00950 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
00951 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
00952 
00953 #ifdef MATLAB_MEX_FILE
00954 #define ASSERT(expression) (mxAssert ((expression), ""))
00955 #else
00956 #define ASSERT(expression) (assert (expression))
00957 #endif /* MATLAB_MEX_FILE */
00958 
00959 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */
00960 (
00961     char *method
00962 ) ;
00963 
00964 PRIVATE void debug_deg_lists
00965 (
00966     Int n_row,
00967     Int n_col,
00968     Colamd_Row Row [],
00969     Colamd_Col Col [],
00970     Int head [],
00971     Int min_score,
00972     Int should,
00973     Int max_deg
00974 ) ;
00975 
00976 PRIVATE void debug_mark
00977 (
00978     Int n_row,
00979     Colamd_Row Row [],
00980     Int tag_mark,
00981     Int max_mark
00982 ) ;
00983 
00984 PRIVATE void debug_matrix
00985 (
00986     Int n_row,
00987     Int n_col,
00988     Colamd_Row Row [],
00989     Colamd_Col Col [],
00990     Int A []
00991 ) ;
00992 
00993 PRIVATE void debug_structures
00994 (
00995     Int n_row,
00996     Int n_col,
00997     Colamd_Row Row [],
00998     Colamd_Col Col [],
00999     Int A [],
01000     Int n_col2
01001 ) ;
01002 
01003 #else /* NDEBUG */
01004 
01005 /* === No debugging ========================================================= */
01006 
01007 #define DEBUG0(params) ;
01008 #define DEBUG1(params) ;
01009 #define DEBUG2(params) ;
01010 #define DEBUG3(params) ;
01011 #define DEBUG4(params) ;
01012 
01013 #define ASSERT(expression)
01014 
01015 #endif /* NDEBUG */
01016 
01017 /* ========================================================================== */
01018 /* === USER-CALLABLE ROUTINES: ============================================== */
01019 /* ========================================================================== */
01020 
01021 /* ========================================================================== */
01022 /* === colamd_recommended =================================================== */
01023 /* ========================================================================== */
01024 
01025 /*
01026     The colamd_recommended routine returns the suggested size for Alen.  This
01027     value has been determined to provide good balance between the number of
01028     garbage collections and the memory requirements for colamd.  If any
01029     argument is negative, or if integer overflow occurs, a 0 is returned as an
01030     error condition.  2*nnz space is required for the row and column
01031     indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
01032     required for the Col and Row arrays, respectively, which are internal to
01033     colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the
01034     minimal amount of "elbow room", and nnz/5 more space is recommended for
01035     run time efficiency.
01036 
01037     Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
01038 
01039     This function is not needed when using symamd.
01040 */
01041 
01042 /* add two values of type size_t, and check for integer overflow */
01043 static size_t t_add (size_t a, size_t b, int *ok)
01044 {
01045     (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
01046     return ((*ok) ? (a + b) : 0) ;
01047 }
01048 
01049 /* compute a*k where k is a small integer, and check for integer overflow */
01050 static size_t t_mult (size_t a, size_t k, int *ok)
01051 {
01052     size_t i, s = 0 ;
01053     for (i = 0 ; i < k ; i++)
01054     {
01055   s = t_add (s, a, ok) ;
01056     }
01057     return (s) ;
01058 }
01059 
01060 /* size of the Col and Row structures */
01061 #define COLAMD_C(n_col,ok) \
01062     ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
01063 
01064 #define COLAMD_R(n_row,ok) \
01065     ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
01066 
01067 
01068 PUBLIC size_t COLAMD_recommended  /* returns recommended value of Alen. */
01069 (
01070     /* === Parameters ======================================================= */
01071 
01072     Int nnz,      /* number of nonzeros in A */
01073     Int n_row,      /* number of rows in A */
01074     Int n_col     /* number of columns in A */
01075 )
01076 {
01077     size_t s, c, r ;
01078     int ok = TRUE ;
01079     if (nnz < 0 || n_row < 0 || n_col < 0)
01080     {
01081   return (0) ;
01082     }
01083     s = t_mult (nnz, 2, &ok) ;      /* 2*nnz */
01084     c = COLAMD_C (n_col, &ok) ;     /* size of column structures */
01085     r = COLAMD_R (n_row, &ok) ;     /* size of row structures */
01086     s = t_add (s, c, &ok) ;
01087     s = t_add (s, r, &ok) ;
01088     s = t_add (s, n_col, &ok) ;     /* elbow room */
01089     s = t_add (s, nnz/5, &ok) ;     /* elbow room */
01090     ok = ok && (s < Int_MAX) ;
01091     return (ok ? s : 0) ;
01092 }
01093 
01094 
01095 /* ========================================================================== */
01096 /* === colamd_set_defaults ================================================== */
01097 /* ========================================================================== */
01098 
01099 /*
01100     The colamd_set_defaults routine sets the default values of the user-
01101     controllable parameters for colamd and symamd:
01102 
01103   Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
01104   entries are removed prior to ordering.  Columns with more than
01105   max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
01106   prior to ordering, and placed last in the output column ordering. 
01107 
01108   Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
01109   entries are removed prior to ordering, and placed last in the
01110   output ordering.
01111 
01112   knobs [0] dense row control
01113 
01114   knobs [1] dense column control
01115 
01116   knobs [2] if nonzero, do aggresive absorption
01117 
01118   knobs [3..19] unused, but future versions might use this
01119 
01120 */
01121 
01122 PUBLIC void COLAMD_set_defaults
01123 (
01124     /* === Parameters ======================================================= */
01125 
01126     double knobs [COLAMD_KNOBS]   /* knob array */
01127 )
01128 {
01129     /* === Local variables ================================================== */
01130 
01131     Int i ;
01132 
01133     if (!knobs)
01134     {
01135   return ;      /* no knobs to initialize */
01136     }
01137     for (i = 0 ; i < COLAMD_KNOBS ; i++)
01138     {
01139   knobs [i] = 0 ;
01140     }
01141     knobs [COLAMD_DENSE_ROW] = 10 ;
01142     knobs [COLAMD_DENSE_COL] = 10 ;
01143     knobs [COLAMD_AGGRESSIVE] = TRUE ;  /* default: do aggressive absorption*/
01144 }
01145 
01146 
01147 /* ========================================================================== */
01148 /* === symamd =============================================================== */
01149 /* ========================================================================== */
01150 
01151 PUBLIC Int SYMAMD_MAIN      /* return TRUE if OK, FALSE otherwise */
01152 (
01153     /* === Parameters ======================================================= */
01154 
01155     Int n,        /* number of rows and columns of A */
01156     Int A [],       /* row indices of A */
01157     Int p [],       /* column pointers of A */
01158     Int perm [],      /* output permutation, size n+1 */
01159     double knobs [COLAMD_KNOBS],  /* parameters (uses defaults if NULL) */
01160     Int stats [COLAMD_STATS],   /* output statistics and error codes */
01161     void * (*allocate) (size_t, size_t),
01162               /* pointer to calloc (ANSI C) or */
01163           /* mxCalloc (for MATLAB mexFunction) */
01164     void (*release) (void *)
01165               /* pointer to free (ANSI C) or */
01166               /* mxFree (for MATLAB mexFunction) */
01167 )
01168 {
01169     /* === Local variables ================================================== */
01170 
01171     Int *count ;    /* length of each column of M, and col pointer*/
01172     Int *mark ;     /* mark array for finding duplicate entries */
01173     Int *M ;      /* row indices of matrix M */
01174     size_t Mlen ;   /* length of M */
01175     Int n_row ;     /* number of rows in M */
01176     Int nnz ;     /* number of entries in A */
01177     Int i ;     /* row index of A */
01178     Int j ;     /* column index of A */
01179     Int k ;     /* row index of M */ 
01180     Int mnz ;     /* number of nonzeros in M */
01181     Int pp ;      /* index into a column of A */
01182     Int last_row ;    /* last row seen in the current column */
01183     Int length ;    /* number of nonzeros in a column */
01184 
01185     double cknobs [COLAMD_KNOBS] ;    /* knobs for colamd */
01186     double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */
01187 
01188 #ifndef NDEBUG
01189     colamd_get_debug ("symamd") ;
01190 #endif /* NDEBUG */
01191 
01192     /* === Check the input arguments ======================================== */
01193 
01194     if (!stats)
01195     {
01196   DEBUG0 (("symamd: stats not present\n")) ;
01197   return (FALSE) ;
01198     }
01199     for (i = 0 ; i < COLAMD_STATS ; i++)
01200     {
01201   stats [i] = 0 ;
01202     }
01203     stats [COLAMD_STATUS] = COLAMD_OK ;
01204     stats [COLAMD_INFO1] = -1 ;
01205     stats [COLAMD_INFO2] = -1 ;
01206 
01207     if (!A)
01208     {
01209       stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
01210   DEBUG0 (("symamd: A not present\n")) ;
01211   return (FALSE) ;
01212     }
01213 
01214     if (!p)   /* p is not present */
01215     {
01216   stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
01217   DEBUG0 (("symamd: p not present\n")) ;
01218       return (FALSE) ;
01219     }
01220 
01221     if (n < 0)    /* n must be >= 0 */
01222     {
01223   stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
01224   stats [COLAMD_INFO1] = n ;
01225   DEBUG0 (("symamd: n negative %d\n", n)) ;
01226       return (FALSE) ;
01227     }
01228 
01229     nnz = p [n] ;
01230     if (nnz < 0)  /* nnz must be >= 0 */
01231     {
01232   stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
01233   stats [COLAMD_INFO1] = nnz ;
01234   DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
01235   return (FALSE) ;
01236     }
01237 
01238     if (p [0] != 0)
01239     {
01240   stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
01241   stats [COLAMD_INFO1] = p [0] ;
01242   DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
01243   return (FALSE) ;
01244     }
01245 
01246     /* === If no knobs, set default knobs =================================== */
01247 
01248     if (!knobs)
01249     {
01250   COLAMD_set_defaults (default_knobs) ;
01251   knobs = default_knobs ;
01252     }
01253 
01254     /* === Allocate count and mark ========================================== */
01255 
01256     count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01257     if (!count)
01258     {
01259   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01260   DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
01261   return (FALSE) ;
01262     }
01263 
01264     mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01265     if (!mark)
01266     {
01267   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01268   (*release) ((void *) count) ;
01269   DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
01270   return (FALSE) ;
01271     }
01272 
01273     /* === Compute column counts of M, check if A is valid ================== */
01274 
01275     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
01276 
01277     for (i = 0 ; i < n ; i++)
01278     {
01279       mark [i] = -1 ;
01280     }
01281 
01282     for (j = 0 ; j < n ; j++)
01283     {
01284   last_row = -1 ;
01285 
01286   length = p [j+1] - p [j] ;
01287   if (length < 0)
01288   {
01289       /* column pointers must be non-decreasing */
01290       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
01291       stats [COLAMD_INFO1] = j ;
01292       stats [COLAMD_INFO2] = length ;
01293       (*release) ((void *) count) ;
01294       (*release) ((void *) mark) ;
01295       DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
01296       return (FALSE) ;
01297   }
01298 
01299   for (pp = p [j] ; pp < p [j+1] ; pp++)
01300   {
01301       i = A [pp] ;
01302       if (i < 0 || i >= n)
01303       {
01304     /* row index i, in column j, is out of bounds */
01305     stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
01306     stats [COLAMD_INFO1] = j ;
01307     stats [COLAMD_INFO2] = i ;
01308     stats [COLAMD_INFO3] = n ;
01309     (*release) ((void *) count) ;
01310     (*release) ((void *) mark) ;
01311     DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
01312     return (FALSE) ;
01313       }
01314 
01315       if (i <= last_row || mark [i] == j)
01316       {
01317     /* row index is unsorted or repeated (or both), thus col */
01318     /* is jumbled.  This is a notice, not an error condition. */
01319     stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
01320     stats [COLAMD_INFO1] = j ;
01321     stats [COLAMD_INFO2] = i ;
01322     (stats [COLAMD_INFO3]) ++ ;
01323     DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
01324       }
01325 
01326       if (i > j && mark [i] != j)
01327       {
01328     /* row k of M will contain column indices i and j */
01329     count [i]++ ;
01330     count [j]++ ;
01331       }
01332 
01333       /* mark the row as having been seen in this column */
01334       mark [i] = j ;
01335 
01336       last_row = i ;
01337   }
01338     }
01339 
01340     /* v2.4: removed free(mark) */
01341 
01342     /* === Compute column pointers of M ===================================== */
01343 
01344     /* use output permutation, perm, for column pointers of M */
01345     perm [0] = 0 ;
01346     for (j = 1 ; j <= n ; j++)
01347     {
01348   perm [j] = perm [j-1] + count [j-1] ;
01349     }
01350     for (j = 0 ; j < n ; j++)
01351     {
01352   count [j] = perm [j] ;
01353     }
01354 
01355     /* === Construct M ====================================================== */
01356 
01357     mnz = perm [n] ;
01358     n_row = mnz / 2 ;
01359     Mlen = COLAMD_recommended (mnz, n_row, n) ;
01360     M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
01361     DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
01362       n_row, n, mnz, (double) Mlen)) ;
01363 
01364     if (!M)
01365     {
01366   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01367   (*release) ((void *) count) ;
01368   (*release) ((void *) mark) ;
01369   DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
01370   return (FALSE) ;
01371     }
01372 
01373     k = 0 ;
01374 
01375     if (stats [COLAMD_STATUS] == COLAMD_OK)
01376     {
01377   /* Matrix is OK */
01378   for (j = 0 ; j < n ; j++)
01379   {
01380       ASSERT (p [j+1] - p [j] >= 0) ;
01381       for (pp = p [j] ; pp < p [j+1] ; pp++)
01382       {
01383     i = A [pp] ;
01384     ASSERT (i >= 0 && i < n) ;
01385     if (i > j)
01386     {
01387         /* row k of M contains column indices i and j */
01388         M [count [i]++] = k ;
01389         M [count [j]++] = k ;
01390         k++ ;
01391     }
01392       }
01393   }
01394     }
01395     else
01396     {
01397   /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
01398   DEBUG0 (("symamd: Duplicates in A.\n")) ;
01399   for (i = 0 ; i < n ; i++)
01400   {
01401       mark [i] = -1 ;
01402   }
01403   for (j = 0 ; j < n ; j++)
01404   {
01405       ASSERT (p [j+1] - p [j] >= 0) ;
01406       for (pp = p [j] ; pp < p [j+1] ; pp++)
01407       {
01408     i = A [pp] ;
01409     ASSERT (i >= 0 && i < n) ;
01410     if (i > j && mark [i] != j)
01411     {
01412         /* row k of M contains column indices i and j */
01413         M [count [i]++] = k ;
01414         M [count [j]++] = k ;
01415         k++ ;
01416         mark [i] = j ;
01417     }
01418       }
01419   }
01420   /* v2.4: free(mark) moved below */
01421     }
01422 
01423     /* count and mark no longer needed */
01424     (*release) ((void *) count) ;
01425     (*release) ((void *) mark) ;  /* v2.4: free (mark) moved here */
01426     ASSERT (k == n_row) ;
01427 
01428     /* === Adjust the knobs for M =========================================== */
01429 
01430     for (i = 0 ; i < COLAMD_KNOBS ; i++)
01431     {
01432   cknobs [i] = knobs [i] ;
01433     }
01434 
01435     /* there are no dense rows in M */
01436     cknobs [COLAMD_DENSE_ROW] = -1 ;
01437     cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
01438 
01439     /* === Order the columns of M =========================================== */
01440 
01441     /* v2.4: colamd cannot fail here, so the error check is removed */
01442     (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
01443 
01444     /* Note that the output permutation is now in perm */
01445 
01446     /* === get the statistics for symamd from colamd ======================== */
01447 
01448     /* a dense column in colamd means a dense row and col in symamd */
01449     stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
01450 
01451     /* === Free M =========================================================== */
01452 
01453     (*release) ((void *) M) ;
01454     DEBUG0 (("symamd: done.\n")) ;
01455     return (TRUE) ;
01456 
01457 }
01458 
01459 /* ========================================================================== */
01460 /* === colamd =============================================================== */
01461 /* ========================================================================== */
01462 
01463 /*
01464     The colamd routine computes a column ordering Q of a sparse matrix
01465     A such that the LU factorization P(AQ) = LU remains sparse, where P is
01466     selected via partial pivoting.   The routine can also be viewed as
01467     providing a permutation Q such that the Cholesky factorization
01468     (AQ)'(AQ) = LL' remains sparse.
01469 */
01470 
01471 PUBLIC Int COLAMD_MAIN    /* returns TRUE if successful, FALSE otherwise*/
01472 (
01473     /* === Parameters ======================================================= */
01474 
01475     Int n_row,      /* number of rows in A */
01476     Int n_col,      /* number of columns in A */
01477     Int Alen,     /* length of A */
01478     Int A [],     /* row indices of A */
01479     Int p [],     /* pointers to columns in A */
01480     double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
01481     Int stats [COLAMD_STATS]  /* output statistics and error codes */
01482 )
01483 {
01484     /* === Local variables ================================================== */
01485 
01486     Int i ;     /* loop index */
01487     Int nnz ;     /* nonzeros in A */
01488     size_t Row_size ;   /* size of Row [], in integers */
01489     size_t Col_size ;   /* size of Col [], in integers */
01490     size_t need ;   /* minimum required length of A */
01491     Colamd_Row *Row ;   /* pointer into A of Row [0..n_row] array */
01492     Colamd_Col *Col ;   /* pointer into A of Col [0..n_col] array */
01493     Int n_col2 ;    /* number of non-dense, non-empty columns */
01494     Int n_row2 ;    /* number of non-dense, non-empty rows */
01495     Int ngarbage ;    /* number of garbage collections performed */
01496     Int max_deg ;   /* maximum row degree */
01497     double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
01498     Int aggressive ;    /* do aggressive absorption */
01499     int ok ;
01500 
01501 #ifndef NDEBUG
01502     colamd_get_debug ("colamd") ;
01503 #endif /* NDEBUG */
01504 
01505     /* === Check the input arguments ======================================== */
01506 
01507     if (!stats)
01508     {
01509   DEBUG0 (("colamd: stats not present\n")) ;
01510   return (FALSE) ;
01511     }
01512     for (i = 0 ; i < COLAMD_STATS ; i++)
01513     {
01514   stats [i] = 0 ;
01515     }
01516     stats [COLAMD_STATUS] = COLAMD_OK ;
01517     stats [COLAMD_INFO1] = -1 ;
01518     stats [COLAMD_INFO2] = -1 ;
01519 
01520     if (!A)   /* A is not present */
01521     {
01522   stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
01523   DEBUG0 (("colamd: A not present\n")) ;
01524   return (FALSE) ;
01525     }
01526 
01527     if (!p)   /* p is not present */
01528     {
01529   stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
01530   DEBUG0 (("colamd: p not present\n")) ;
01531       return (FALSE) ;
01532     }
01533 
01534     if (n_row < 0)  /* n_row must be >= 0 */
01535     {
01536   stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
01537   stats [COLAMD_INFO1] = n_row ;
01538   DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
01539       return (FALSE) ;
01540     }
01541 
01542     if (n_col < 0)  /* n_col must be >= 0 */
01543     {
01544   stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
01545   stats [COLAMD_INFO1] = n_col ;
01546   DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
01547       return (FALSE) ;
01548     }
01549 
01550     nnz = p [n_col] ;
01551     if (nnz < 0)  /* nnz must be >= 0 */
01552     {
01553   stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
01554   stats [COLAMD_INFO1] = nnz ;
01555   DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
01556   return (FALSE) ;
01557     }
01558 
01559     if (p [0] != 0)
01560     {
01561   stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
01562   stats [COLAMD_INFO1] = p [0] ;
01563   DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
01564   return (FALSE) ;
01565     }
01566 
01567     /* === If no knobs, set default knobs =================================== */
01568 
01569     if (!knobs)
01570     {
01571   COLAMD_set_defaults (default_knobs) ;
01572   knobs = default_knobs ;
01573     }
01574 
01575     aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
01576 
01577     /* === Allocate the Row and Col arrays from array A ===================== */
01578 
01579     ok = TRUE ;
01580     Col_size = COLAMD_C (n_col, &ok) ;      /* size of Col array of structs */
01581     Row_size = COLAMD_R (n_row, &ok) ;      /* size of Row array of structs */
01582 
01583     /* need = 2*nnz + n_col + Col_size + Row_size ; */
01584     need = t_mult (nnz, 2, &ok) ;
01585     need = t_add (need, n_col, &ok) ;
01586     need = t_add (need, Col_size, &ok) ;
01587     need = t_add (need, Row_size, &ok) ;
01588 
01589     if (!ok || need > (size_t) Alen || need > Int_MAX)
01590     {
01591   /* not enough space in array A to perform the ordering */
01592   stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
01593   stats [COLAMD_INFO1] = need ;
01594   stats [COLAMD_INFO2] = Alen ;
01595   DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
01596   return (FALSE) ;
01597     }
01598 
01599     Alen -= Col_size + Row_size ;
01600     Col = (Colamd_Col *) &A [Alen] ;
01601     Row = (Colamd_Row *) &A [Alen + Col_size] ;
01602 
01603     /* === Construct the row and column data structures ===================== */
01604 
01605     if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
01606     {
01607   /* input matrix is invalid */
01608   DEBUG0 (("colamd: Matrix invalid\n")) ;
01609   return (FALSE) ;
01610     }
01611 
01612     /* === Initialize scores, kill dense rows/columns ======================= */
01613 
01614     init_scoring (n_row, n_col, Row, Col, A, p, knobs,
01615   &n_row2, &n_col2, &max_deg) ;
01616 
01617     /* === Order the supercolumns =========================================== */
01618 
01619     ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
01620   n_col2, max_deg, 2*nnz, aggressive) ;
01621 
01622     /* === Order the non-principal columns ================================== */
01623 
01624     order_children (n_col, Col, p) ;
01625 
01626     /* === Return statistics in stats ======================================= */
01627 
01628     stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
01629     stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
01630     stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
01631     DEBUG0 (("colamd: done.\n")) ; 
01632     return (TRUE) ;
01633 }
01634 
01635 
01636 /* ========================================================================== */
01637 /* === colamd_report ======================================================== */
01638 /* ========================================================================== */
01639 
01640 PUBLIC void COLAMD_report
01641 (
01642     Int stats [COLAMD_STATS]
01643 )
01644 {
01645     print_report ("colamd", stats) ;
01646 }
01647 
01648 
01649 /* ========================================================================== */
01650 /* === symamd_report ======================================================== */
01651 /* ========================================================================== */
01652 
01653 PUBLIC void SYMAMD_report
01654 (
01655     Int stats [COLAMD_STATS]
01656 )
01657 {
01658     print_report ("symamd", stats) ;
01659 }
01660 
01661 
01662 
01663 /* ========================================================================== */
01664 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
01665 /* ========================================================================== */
01666 
01667 /* There are no user-callable routines beyond this point in the file */
01668 
01669 
01670 /* ========================================================================== */
01671 /* === init_rows_cols ======================================================= */
01672 /* ========================================================================== */
01673 
01674 /*
01675     Takes the column form of the matrix in A and creates the row form of the
01676     matrix.  Also, row and column attributes are stored in the Col and Row
01677     structs.  If the columns are un-sorted or contain duplicate row indices,
01678     this routine will also sort and remove duplicate row indices from the
01679     column form of the matrix.  Returns FALSE if the matrix is invalid,
01680     TRUE otherwise.  Not user-callable.
01681 */
01682 
01683 PRIVATE Int init_rows_cols  /* returns TRUE if OK, or FALSE otherwise */
01684 (
01685     /* === Parameters ======================================================= */
01686 
01687     Int n_row,      /* number of rows of A */
01688     Int n_col,      /* number of columns of A */
01689     Colamd_Row Row [],    /* of size n_row+1 */
01690     Colamd_Col Col [],    /* of size n_col+1 */
01691     Int A [],     /* row indices of A, of size Alen */
01692     Int p [],     /* pointers to columns in A, of size n_col+1 */
01693     Int stats [COLAMD_STATS]  /* colamd statistics */ 
01694 )
01695 {
01696     /* === Local variables ================================================== */
01697 
01698     Int col ;     /* a column index */
01699     Int row ;     /* a row index */
01700     Int *cp ;     /* a column pointer */
01701     Int *cp_end ;   /* a pointer to the end of a column */
01702     Int *rp ;     /* a row pointer */
01703     Int *rp_end ;   /* a pointer to the end of a row */
01704     Int last_row ;    /* previous row */
01705 
01706     /* === Initialize columns, and check column pointers ==================== */
01707 
01708     for (col = 0 ; col < n_col ; col++)
01709     {
01710   Col [col].start = p [col] ;
01711   Col [col].length = p [col+1] - p [col] ;
01712 
01713   if (Col [col].length < 0)
01714   {
01715       /* column pointers must be non-decreasing */
01716       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
01717       stats [COLAMD_INFO1] = col ;
01718       stats [COLAMD_INFO2] = Col [col].length ;
01719       DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
01720       return (FALSE) ;
01721   }
01722 
01723   Col [col].shared1.thickness = 1 ;
01724   Col [col].shared2.score = 0 ;
01725   Col [col].shared3.prev = EMPTY ;
01726   Col [col].shared4.degree_next = EMPTY ;
01727     }
01728 
01729     /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
01730 
01731     /* === Scan columns, compute row degrees, and check row indices ========= */
01732 
01733     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
01734 
01735     for (row = 0 ; row < n_row ; row++)
01736     {
01737   Row [row].length = 0 ;
01738   Row [row].shared2.mark = -1 ;
01739     }
01740 
01741     for (col = 0 ; col < n_col ; col++)
01742     {
01743   last_row = -1 ;
01744 
01745   cp = &A [p [col]] ;
01746   cp_end = &A [p [col+1]] ;
01747 
01748   while (cp < cp_end)
01749   {
01750       row = *cp++ ;
01751 
01752       /* make sure row indices within range */
01753       if (row < 0 || row >= n_row)
01754       {
01755     stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
01756     stats [COLAMD_INFO1] = col ;
01757     stats [COLAMD_INFO2] = row ;
01758     stats [COLAMD_INFO3] = n_row ;
01759     DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
01760     return (FALSE) ;
01761       }
01762 
01763       if (row <= last_row || Row [row].shared2.mark == col)
01764       {
01765     /* row index are unsorted or repeated (or both), thus col */
01766     /* is jumbled.  This is a notice, not an error condition. */
01767     stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
01768     stats [COLAMD_INFO1] = col ;
01769     stats [COLAMD_INFO2] = row ;
01770     (stats [COLAMD_INFO3]) ++ ;
01771     DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
01772       }
01773 
01774       if (Row [row].shared2.mark != col)
01775       {
01776     Row [row].length++ ;
01777       }
01778       else
01779       {
01780     /* this is a repeated entry in the column, */
01781     /* it will be removed */
01782     Col [col].length-- ;
01783       }
01784 
01785       /* mark the row as having been seen in this column */
01786       Row [row].shared2.mark = col ;
01787 
01788       last_row = row ;
01789   }
01790     }
01791 
01792     /* === Compute row pointers ============================================= */
01793 
01794     /* row form of the matrix starts directly after the column */
01795     /* form of matrix in A */
01796     Row [0].start = p [n_col] ;
01797     Row [0].shared1.p = Row [0].start ;
01798     Row [0].shared2.mark = -1 ;
01799     for (row = 1 ; row < n_row ; row++)
01800     {
01801   Row [row].start = Row [row-1].start + Row [row-1].length ;
01802   Row [row].shared1.p = Row [row].start ;
01803   Row [row].shared2.mark = -1 ;
01804     }
01805 
01806     /* === Create row form ================================================== */
01807 
01808     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
01809     {
01810   /* if cols jumbled, watch for repeated row indices */
01811   for (col = 0 ; col < n_col ; col++)
01812   {
01813       cp = &A [p [col]] ;
01814       cp_end = &A [p [col+1]] ;
01815       while (cp < cp_end)
01816       {
01817     row = *cp++ ;
01818     if (Row [row].shared2.mark != col)
01819     {
01820         A [(Row [row].shared1.p)++] = col ;
01821         Row [row].shared2.mark = col ;
01822     }
01823       }
01824   }
01825     }
01826     else
01827     {
01828   /* if cols not jumbled, we don't need the mark (this is faster) */
01829   for (col = 0 ; col < n_col ; col++)
01830   {
01831       cp = &A [p [col]] ;
01832       cp_end = &A [p [col+1]] ;
01833       while (cp < cp_end)
01834       {
01835     A [(Row [*cp++].shared1.p)++] = col ;
01836       }
01837   }
01838     }
01839 
01840     /* === Clear the row marks and set row degrees ========================== */
01841 
01842     for (row = 0 ; row < n_row ; row++)
01843     {
01844   Row [row].shared2.mark = 0 ;
01845   Row [row].shared1.degree = Row [row].length ;
01846     }
01847 
01848     /* === See if we need to re-create columns ============================== */
01849 
01850     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
01851     {
01852       DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
01853 
01854 #ifndef NDEBUG
01855   /* make sure column lengths are correct */
01856   for (col = 0 ; col < n_col ; col++)
01857   {
01858       p [col] = Col [col].length ;
01859   }
01860   for (row = 0 ; row < n_row ; row++)
01861   {
01862       rp = &A [Row [row].start] ;
01863       rp_end = rp + Row [row].length ;
01864       while (rp < rp_end)
01865       {
01866     p [*rp++]-- ;
01867       }
01868   }
01869   for (col = 0 ; col < n_col ; col++)
01870   {
01871       ASSERT (p [col] == 0) ;
01872   }
01873   /* now p is all zero (different than when debugging is turned off) */
01874 #endif /* NDEBUG */
01875 
01876   /* === Compute col pointers ========================================= */
01877 
01878   /* col form of the matrix starts at A [0]. */
01879   /* Note, we may have a gap between the col form and the row */
01880   /* form if there were duplicate entries, if so, it will be */
01881   /* removed upon the first garbage collection */
01882   Col [0].start = 0 ;
01883   p [0] = Col [0].start ;
01884   for (col = 1 ; col < n_col ; col++)
01885   {
01886       /* note that the lengths here are for pruned columns, i.e. */
01887       /* no duplicate row indices will exist for these columns */
01888       Col [col].start = Col [col-1].start + Col [col-1].length ;
01889       p [col] = Col [col].start ;
01890   }
01891 
01892   /* === Re-create col form =========================================== */
01893 
01894   for (row = 0 ; row < n_row ; row++)
01895   {
01896       rp = &A [Row [row].start] ;
01897       rp_end = rp + Row [row].length ;
01898       while (rp < rp_end)
01899       {
01900     A [(p [*rp++])++] = row ;
01901       }
01902   }
01903     }
01904 
01905     /* === Done.  Matrix is not (or no longer) jumbled ====================== */
01906 
01907     return (TRUE) ;
01908 }
01909 
01910 
01911 /* ========================================================================== */
01912 /* === init_scoring ========================================================= */
01913 /* ========================================================================== */
01914 
01915 /*
01916     Kills dense or empty columns and rows, calculates an initial score for
01917     each column, and places all columns in the degree lists.  Not user-callable.
01918 */
01919 
01920 PRIVATE void init_scoring
01921 (
01922     /* === Parameters ======================================================= */
01923 
01924     Int n_row,      /* number of rows of A */
01925     Int n_col,      /* number of columns of A */
01926     Colamd_Row Row [],    /* of size n_row+1 */
01927     Colamd_Col Col [],    /* of size n_col+1 */
01928     Int A [],     /* column form and row form of A */
01929     Int head [],    /* of size n_col+1 */
01930     double knobs [COLAMD_KNOBS],/* parameters */
01931     Int *p_n_row2,    /* number of non-dense, non-empty rows */
01932     Int *p_n_col2,    /* number of non-dense, non-empty columns */
01933     Int *p_max_deg    /* maximum row degree */
01934 )
01935 {
01936     /* === Local variables ================================================== */
01937 
01938     Int c ;     /* a column index */
01939     Int r, row ;    /* a row index */
01940     Int *cp ;     /* a column pointer */
01941     Int deg ;     /* degree of a row or column */
01942     Int *cp_end ;   /* a pointer to the end of a column */
01943     Int *new_cp ;   /* new column pointer */
01944     Int col_length ;    /* length of pruned column */
01945     Int score ;     /* current column score */
01946     Int n_col2 ;    /* number of non-dense, non-empty columns */
01947     Int n_row2 ;    /* number of non-dense, non-empty rows */
01948     Int dense_row_count ; /* remove rows with more entries than this */
01949     Int dense_col_count ; /* remove cols with more entries than this */
01950     Int min_score ;   /* smallest column score */
01951     Int max_deg ;   /* maximum row degree */
01952     Int next_col ;    /* Used to add to degree list.*/
01953 
01954 #ifndef NDEBUG
01955     Int debug_count ;   /* debug only. */
01956 #endif /* NDEBUG */
01957 
01958     /* === Extract knobs ==================================================== */
01959 
01960     /* Note: if knobs contains a NaN, this is undefined: */
01961     if (knobs [COLAMD_DENSE_ROW] < 0)
01962     {
01963   /* only remove completely dense rows */
01964   dense_row_count = n_col-1 ;
01965     }
01966     else
01967     {
01968   dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
01969     }
01970     if (knobs [COLAMD_DENSE_COL] < 0)
01971     {
01972   /* only remove completely dense columns */
01973   dense_col_count = n_row-1 ;
01974     }
01975     else
01976     {
01977   dense_col_count =
01978       DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
01979     }
01980 
01981     DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
01982     max_deg = 0 ;
01983     n_col2 = n_col ;
01984     n_row2 = n_row ;
01985 
01986     /* === Kill empty columns =============================================== */
01987 
01988     /* Put the empty columns at the end in their natural order, so that LU */
01989     /* factorization can proceed as far as possible. */
01990     for (c = n_col-1 ; c >= 0 ; c--)
01991     {
01992   deg = Col [c].length ;
01993   if (deg == 0)
01994   {
01995       /* this is a empty column, kill and order it last */
01996       Col [c].shared2.order = --n_col2 ;
01997       KILL_PRINCIPAL_COL (c) ;
01998   }
01999     }
02000     DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
02001 
02002     /* === Kill dense columns =============================================== */
02003 
02004     /* Put the dense columns at the end, in their natural order */
02005     for (c = n_col-1 ; c >= 0 ; c--)
02006     {
02007   /* skip any dead columns */
02008   if (COL_IS_DEAD (c))
02009   {
02010       continue ;
02011   }
02012   deg = Col [c].length ;
02013   if (deg > dense_col_count)
02014   {
02015       /* this is a dense column, kill and order it last */
02016       Col [c].shared2.order = --n_col2 ;
02017       /* decrement the row degrees */
02018       cp = &A [Col [c].start] ;
02019       cp_end = cp + Col [c].length ;
02020       while (cp < cp_end)
02021       {
02022     Row [*cp++].shared1.degree-- ;
02023       }
02024       KILL_PRINCIPAL_COL (c) ;
02025   }
02026     }
02027     DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
02028 
02029     /* === Kill dense and empty rows ======================================== */
02030 
02031     for (r = 0 ; r < n_row ; r++)
02032     {
02033   deg = Row [r].shared1.degree ;
02034   ASSERT (deg >= 0 && deg <= n_col) ;
02035   if (deg > dense_row_count || deg == 0)
02036   {
02037       /* kill a dense or empty row */
02038       KILL_ROW (r) ;
02039       --n_row2 ;
02040   }
02041   else
02042   {
02043       /* keep track of max degree of remaining rows */
02044       max_deg = MAX (max_deg, deg) ;
02045   }
02046     }
02047     DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
02048 
02049     /* === Compute initial column scores ==================================== */
02050 
02051     /* At this point the row degrees are accurate.  They reflect the number */
02052     /* of "live" (non-dense) columns in each row.  No empty rows exist. */
02053     /* Some "live" columns may contain only dead rows, however.  These are */
02054     /* pruned in the code below. */
02055 
02056     /* now find the initial matlab score for each column */
02057     for (c = n_col-1 ; c >= 0 ; c--)
02058     {
02059   /* skip dead column */
02060   if (COL_IS_DEAD (c))
02061   {
02062       continue ;
02063   }
02064   score = 0 ;
02065   cp = &A [Col [c].start] ;
02066   new_cp = cp ;
02067   cp_end = cp + Col [c].length ;
02068   while (cp < cp_end)
02069   {
02070       /* get a row */
02071       row = *cp++ ;
02072       /* skip if dead */
02073       if (ROW_IS_DEAD (row))
02074       {
02075     continue ;
02076       }
02077       /* compact the column */
02078       *new_cp++ = row ;
02079       /* add row's external degree */
02080       score += Row [row].shared1.degree - 1 ;
02081       /* guard against integer overflow */
02082       score = MIN (score, n_col) ;
02083   }
02084   /* determine pruned column length */
02085   col_length = (Int) (new_cp - &A [Col [c].start]) ;
02086   if (col_length == 0)
02087   {
02088       /* a newly-made null column (all rows in this col are "dense" */
02089       /* and have already been killed) */
02090       DEBUG2 (("Newly null killed: %d\n", c)) ;
02091       Col [c].shared2.order = --n_col2 ;
02092       KILL_PRINCIPAL_COL (c) ;
02093   }
02094   else
02095   {
02096       /* set column length and set score */
02097       ASSERT (score >= 0) ;
02098       ASSERT (score <= n_col) ;
02099       Col [c].length = col_length ;
02100       Col [c].shared2.score = score ;
02101   }
02102     }
02103     DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
02104       n_col-n_col2)) ;
02105 
02106     /* At this point, all empty rows and columns are dead.  All live columns */
02107     /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
02108     /* yet).  Rows may contain dead columns, but all live rows contain at */
02109     /* least one live column. */
02110 
02111 #ifndef NDEBUG
02112     debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
02113 #endif /* NDEBUG */
02114 
02115     /* === Initialize degree lists ========================================== */
02116 
02117 #ifndef NDEBUG
02118     debug_count = 0 ;
02119 #endif /* NDEBUG */
02120 
02121     /* clear the hash buckets */
02122     for (c = 0 ; c <= n_col ; c++)
02123     {
02124   head [c] = EMPTY ;
02125     }
02126     min_score = n_col ;
02127     /* place in reverse order, so low column indices are at the front */
02128     /* of the lists.  This is to encourage natural tie-breaking */
02129     for (c = n_col-1 ; c >= 0 ; c--)
02130     {
02131   /* only add principal columns to degree lists */
02132   if (COL_IS_ALIVE (c))
02133   {
02134       DEBUG4 (("place %d score %d minscore %d ncol %d\n",
02135     c, Col [c].shared2.score, min_score, n_col)) ;
02136 
02137       /* === Add columns score to DList =============================== */
02138 
02139       score = Col [c].shared2.score ;
02140 
02141       ASSERT (min_score >= 0) ;
02142       ASSERT (min_score <= n_col) ;
02143       ASSERT (score >= 0) ;
02144       ASSERT (score <= n_col) ;
02145       ASSERT (head [score] >= EMPTY) ;
02146 
02147       /* now add this column to dList at proper score location */
02148       next_col = head [score] ;
02149       Col [c].shared3.prev = EMPTY ;
02150       Col [c].shared4.degree_next = next_col ;
02151 
02152       /* if there already was a column with the same score, set its */
02153       /* previous pointer to this new column */
02154       if (next_col != EMPTY)
02155       {
02156     Col [next_col].shared3.prev = c ;
02157       }
02158       head [score] = c ;
02159 
02160       /* see if this score is less than current min */
02161       min_score = MIN (min_score, score) ;
02162 
02163 #ifndef NDEBUG
02164       debug_count++ ;
02165 #endif /* NDEBUG */
02166 
02167   }
02168     }
02169 
02170 #ifndef NDEBUG
02171     DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
02172   debug_count, n_col, n_col-debug_count)) ;
02173     ASSERT (debug_count == n_col2) ;
02174     debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
02175 #endif /* NDEBUG */
02176 
02177     /* === Return number of remaining columns, and max row degree =========== */
02178 
02179     *p_n_col2 = n_col2 ;
02180     *p_n_row2 = n_row2 ;
02181     *p_max_deg = max_deg ;
02182 }
02183 
02184 
02185 /* ========================================================================== */
02186 /* === find_ordering ======================================================== */
02187 /* ========================================================================== */
02188 
02189 /*
02190     Order the principal columns of the supercolumn form of the matrix
02191     (no supercolumns on input).  Uses a minimum approximate column minimum
02192     degree ordering method.  Not user-callable.
02193 */
02194 
02195 PRIVATE Int find_ordering /* return the number of garbage collections */
02196 (
02197     /* === Parameters ======================================================= */
02198 
02199     Int n_row,      /* number of rows of A */
02200     Int n_col,      /* number of columns of A */
02201     Int Alen,     /* size of A, 2*nnz + n_col or larger */
02202     Colamd_Row Row [],    /* of size n_row+1 */
02203     Colamd_Col Col [],    /* of size n_col+1 */
02204     Int A [],     /* column form and row form of A */
02205     Int head [],    /* of size n_col+1 */
02206     Int n_col2,     /* Remaining columns to order */
02207     Int max_deg,    /* Maximum row degree */
02208     Int pfree,      /* index of first free slot (2*nnz on entry) */
02209     Int aggressive
02210 )
02211 {
02212     /* === Local variables ================================================== */
02213 
02214     Int k ;     /* current pivot ordering step */
02215     Int pivot_col ;   /* current pivot column */
02216     Int *cp ;     /* a column pointer */
02217     Int *rp ;     /* a row pointer */
02218     Int pivot_row ;   /* current pivot row */
02219     Int *new_cp ;   /* modified column pointer */
02220     Int *new_rp ;   /* modified row pointer */
02221     Int pivot_row_start ; /* pointer to start of pivot row */
02222     Int pivot_row_degree ;  /* number of columns in pivot row */
02223     Int pivot_row_length ;  /* number of supercolumns in pivot row */
02224     Int pivot_col_score ; /* score of pivot column */
02225     Int needed_memory ;   /* free space needed for pivot row */
02226     Int *cp_end ;   /* pointer to the end of a column */
02227     Int *rp_end ;   /* pointer to the end of a row */
02228     Int row ;     /* a row index */
02229     Int col ;     /* a column index */
02230     Int max_score ;   /* maximum possible score */
02231     Int cur_score ;   /* score of current column */
02232     unsigned Int hash ;   /* hash value for supernode detection */
02233     Int head_column ;   /* head of hash bucket */
02234     Int first_col ;   /* first column in hash bucket */
02235     Int tag_mark ;    /* marker value for mark array */
02236     Int row_mark ;    /* Row [row].shared2.mark */
02237     Int set_difference ;  /* set difference size of row with pivot row */
02238     Int min_score ;   /* smallest column score */
02239     Int col_thickness ;   /* "thickness" (no. of columns in a supercol) */
02240     Int max_mark ;    /* maximum value of tag_mark */
02241     Int pivot_col_thickness ; /* number of columns represented by pivot col */
02242     Int prev_col ;    /* Used by Dlist operations. */
02243     Int next_col ;    /* Used by Dlist operations. */
02244     Int ngarbage ;    /* number of garbage collections performed */
02245 
02246 #ifndef NDEBUG
02247     Int debug_d ;   /* debug loop counter */
02248     Int debug_step = 0 ;  /* debug loop counter */
02249 #endif /* NDEBUG */
02250 
02251     /* === Initialization and clear mark ==================================== */
02252 
02253     max_mark = INT_MAX - n_col ;  /* INT_MAX defined in <limits.h> */
02254     tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02255     min_score = 0 ;
02256     ngarbage = 0 ;
02257     DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
02258 
02259     /* === Order the columns ================================================ */
02260 
02261     for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
02262     {
02263 
02264 #ifndef NDEBUG
02265   if (debug_step % 100 == 0)
02266   {
02267       DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
02268   }
02269   else
02270   {
02271       DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
02272   }
02273   debug_step++ ;
02274   debug_deg_lists (n_row, n_col, Row, Col, head,
02275     min_score, n_col2-k, max_deg) ;
02276   debug_matrix (n_row, n_col, Row, Col, A) ;
02277 #endif /* NDEBUG */
02278 
02279   /* === Select pivot column, and order it ============================ */
02280 
02281   /* make sure degree list isn't empty */
02282   ASSERT (min_score >= 0) ;
02283   ASSERT (min_score <= n_col) ;
02284   ASSERT (head [min_score] >= EMPTY) ;
02285 
02286 #ifndef NDEBUG
02287   for (debug_d = 0 ; debug_d < min_score ; debug_d++)
02288   {
02289       ASSERT (head [debug_d] == EMPTY) ;
02290   }
02291 #endif /* NDEBUG */
02292 
02293   /* get pivot column from head of minimum degree list */
02294   while (head [min_score] == EMPTY && min_score < n_col)
02295   {
02296       min_score++ ;
02297   }
02298   pivot_col = head [min_score] ;
02299   ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
02300   next_col = Col [pivot_col].shared4.degree_next ;
02301   head [min_score] = next_col ;
02302   if (next_col != EMPTY)
02303   {
02304       Col [next_col].shared3.prev = EMPTY ;
02305   }
02306 
02307   ASSERT (COL_IS_ALIVE (pivot_col)) ;
02308 
02309   /* remember score for defrag check */
02310   pivot_col_score = Col [pivot_col].shared2.score ;
02311 
02312   /* the pivot column is the kth column in the pivot order */
02313   Col [pivot_col].shared2.order = k ;
02314 
02315   /* increment order count by column thickness */
02316   pivot_col_thickness = Col [pivot_col].shared1.thickness ;
02317   k += pivot_col_thickness ;
02318   ASSERT (pivot_col_thickness > 0) ;
02319   DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
02320 
02321   /* === Garbage_collection, if necessary ============================= */
02322 
02323   needed_memory = MIN (pivot_col_score, n_col - k) ;
02324   if (pfree + needed_memory >= Alen)
02325   {
02326       pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
02327       ngarbage++ ;
02328       /* after garbage collection we will have enough */
02329       ASSERT (pfree + needed_memory < Alen) ;
02330       /* garbage collection has wiped out the Row[].shared2.mark array */
02331       tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02332 
02333 #ifndef NDEBUG
02334       debug_matrix (n_row, n_col, Row, Col, A) ;
02335 #endif /* NDEBUG */
02336   }
02337 
02338   /* === Compute pivot row pattern ==================================== */
02339 
02340   /* get starting location for this new merged row */
02341   pivot_row_start = pfree ;
02342 
02343   /* initialize new row counts to zero */
02344   pivot_row_degree = 0 ;
02345 
02346   /* tag pivot column as having been visited so it isn't included */
02347   /* in merged pivot row */
02348   Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
02349 
02350   /* pivot row is the union of all rows in the pivot column pattern */
02351   cp = &A [Col [pivot_col].start] ;
02352   cp_end = cp + Col [pivot_col].length ;
02353   while (cp < cp_end)
02354   {
02355       /* get a row */
02356       row = *cp++ ;
02357       DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
02358       /* skip if row is dead */
02359       if (ROW_IS_ALIVE (row))
02360       {
02361     rp = &A [Row [row].start] ;
02362     rp_end = rp + Row [row].length ;
02363     while (rp < rp_end)
02364     {
02365         /* get a column */
02366         col = *rp++ ;
02367         /* add the column, if alive and untagged */
02368         col_thickness = Col [col].shared1.thickness ;
02369         if (col_thickness > 0 && COL_IS_ALIVE (col))
02370         {
02371       /* tag column in pivot row */
02372       Col [col].shared1.thickness = -col_thickness ;
02373       ASSERT (pfree < Alen) ;
02374       /* place column in pivot row */
02375       A [pfree++] = col ;
02376       pivot_row_degree += col_thickness ;
02377         }
02378     }
02379       }
02380   }
02381 
02382   /* clear tag on pivot column */
02383   Col [pivot_col].shared1.thickness = pivot_col_thickness ;
02384   max_deg = MAX (max_deg, pivot_row_degree) ;
02385 
02386 #ifndef NDEBUG
02387   DEBUG3 (("check2\n")) ;
02388   debug_mark (n_row, Row, tag_mark, max_mark) ;
02389 #endif /* NDEBUG */
02390 
02391   /* === Kill all rows used to construct pivot row ==================== */
02392 
02393   /* also kill pivot row, temporarily */
02394   cp = &A [Col [pivot_col].start] ;
02395   cp_end = cp + Col [pivot_col].length ;
02396   while (cp < cp_end)
02397   {
02398       /* may be killing an already dead row */
02399       row = *cp++ ;
02400       DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
02401       KILL_ROW (row) ;
02402   }
02403 
02404   /* === Select a row index to use as the new pivot row =============== */
02405 
02406   pivot_row_length = pfree - pivot_row_start ;
02407   if (pivot_row_length > 0)
02408   {
02409       /* pick the "pivot" row arbitrarily (first row in col) */
02410       pivot_row = A [Col [pivot_col].start] ;
02411       DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
02412   }
02413   else
02414   {
02415       /* there is no pivot row, since it is of zero length */
02416       pivot_row = EMPTY ;
02417       ASSERT (pivot_row_length == 0) ;
02418   }
02419   ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
02420 
02421   /* === Approximate degree computation =============================== */
02422 
02423   /* Here begins the computation of the approximate degree.  The column */
02424   /* score is the sum of the pivot row "length", plus the size of the */
02425   /* set differences of each row in the column minus the pattern of the */
02426   /* pivot row itself.  The column ("thickness") itself is also */
02427   /* excluded from the column score (we thus use an approximate */
02428   /* external degree). */
02429 
02430   /* The time taken by the following code (compute set differences, and */
02431   /* add them up) is proportional to the size of the data structure */
02432   /* being scanned - that is, the sum of the sizes of each column in */
02433   /* the pivot row.  Thus, the amortized time to compute a column score */
02434   /* is proportional to the size of that column (where size, in this */
02435   /* context, is the column "length", or the number of row indices */
02436   /* in that column).  The number of row indices in a column is */
02437   /* monotonically non-decreasing, from the length of the original */
02438   /* column on input to colamd. */
02439 
02440   /* === Compute set differences ====================================== */
02441 
02442   DEBUG3 (("** Computing set differences phase. **\n")) ;
02443 
02444   /* pivot row is currently dead - it will be revived later. */
02445 
02446   DEBUG3 (("Pivot row: ")) ;
02447   /* for each column in pivot row */
02448   rp = &A [pivot_row_start] ;
02449   rp_end = rp + pivot_row_length ;
02450   while (rp < rp_end)
02451   {
02452       col = *rp++ ;
02453       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
02454       DEBUG3 (("Col: %d\n", col)) ;
02455 
02456       /* clear tags used to construct pivot row pattern */
02457       col_thickness = -Col [col].shared1.thickness ;
02458       ASSERT (col_thickness > 0) ;
02459       Col [col].shared1.thickness = col_thickness ;
02460 
02461       /* === Remove column from degree list =========================== */
02462 
02463       cur_score = Col [col].shared2.score ;
02464       prev_col = Col [col].shared3.prev ;
02465       next_col = Col [col].shared4.degree_next ;
02466       ASSERT (cur_score >= 0) ;
02467       ASSERT (cur_score <= n_col) ;
02468       ASSERT (cur_score >= EMPTY) ;
02469       if (prev_col == EMPTY)
02470       {
02471     head [cur_score] = next_col ;
02472       }
02473       else
02474       {
02475     Col [prev_col].shared4.degree_next = next_col ;
02476       }
02477       if (next_col != EMPTY)
02478       {
02479     Col [next_col].shared3.prev = prev_col ;
02480       }
02481 
02482       /* === Scan the column ========================================== */
02483 
02484       cp = &A [Col [col].start] ;
02485       cp_end = cp + Col [col].length ;
02486       while (cp < cp_end)
02487       {
02488     /* get a row */
02489     row = *cp++ ;
02490     row_mark = Row [row].shared2.mark ;
02491     /* skip if dead */
02492     if (ROW_IS_MARKED_DEAD (row_mark))
02493     {
02494         continue ;
02495     }
02496     ASSERT (row != pivot_row) ;
02497     set_difference = row_mark - tag_mark ;
02498     /* check if the row has been seen yet */
02499     if (set_difference < 0)
02500     {
02501         ASSERT (Row [row].shared1.degree <= max_deg) ;
02502         set_difference = Row [row].shared1.degree ;
02503     }
02504     /* subtract column thickness from this row's set difference */
02505     set_difference -= col_thickness ;
02506     ASSERT (set_difference >= 0) ;
02507     /* absorb this row if the set difference becomes zero */
02508     if (set_difference == 0 && aggressive)
02509     {
02510         DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
02511         KILL_ROW (row) ;
02512     }
02513     else
02514     {
02515         /* save the new mark */
02516         Row [row].shared2.mark = set_difference + tag_mark ;
02517     }
02518       }
02519   }
02520 
02521 #ifndef NDEBUG
02522   debug_deg_lists (n_row, n_col, Row, Col, head,
02523     min_score, n_col2-k-pivot_row_degree, max_deg) ;
02524 #endif /* NDEBUG */
02525 
02526   /* === Add up set differences for each column ======================= */
02527 
02528   DEBUG3 (("** Adding set differences phase. **\n")) ;
02529 
02530   /* for each column in pivot row */
02531   rp = &A [pivot_row_start] ;
02532   rp_end = rp + pivot_row_length ;
02533   while (rp < rp_end)
02534   {
02535       /* get a column */
02536       col = *rp++ ;
02537       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
02538       hash = 0 ;
02539       cur_score = 0 ;
02540       cp = &A [Col [col].start] ;
02541       /* compact the column */
02542       new_cp = cp ;
02543       cp_end = cp + Col [col].length ;
02544 
02545       DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
02546 
02547       while (cp < cp_end)
02548       {
02549     /* get a row */
02550     row = *cp++ ;
02551     ASSERT(row >= 0 && row < n_row) ;
02552     row_mark = Row [row].shared2.mark ;
02553     /* skip if dead */
02554     if (ROW_IS_MARKED_DEAD (row_mark))
02555     {
02556         DEBUG4 ((" Row %d, dead\n", row)) ;
02557         continue ;
02558     }
02559     DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
02560     ASSERT (row_mark >= tag_mark) ;
02561     /* compact the column */
02562     *new_cp++ = row ;
02563     /* compute hash function */
02564     hash += row ;
02565     /* add set difference */
02566     cur_score += row_mark - tag_mark ;
02567     /* integer overflow... */
02568     cur_score = MIN (cur_score, n_col) ;
02569       }
02570 
02571       /* recompute the column's length */
02572       Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
02573 
02574       /* === Further mass elimination ================================= */
02575 
02576       if (Col [col].length == 0)
02577       {
02578     DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
02579     /* nothing left but the pivot row in this column */
02580     KILL_PRINCIPAL_COL (col) ;
02581     pivot_row_degree -= Col [col].shared1.thickness ;
02582     ASSERT (pivot_row_degree >= 0) ;
02583     /* order it */
02584     Col [col].shared2.order = k ;
02585     /* increment order count by column thickness */
02586     k += Col [col].shared1.thickness ;
02587       }
02588       else
02589       {
02590     /* === Prepare for supercolumn detection ==================== */
02591 
02592     DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
02593 
02594     /* save score so far */
02595     Col [col].shared2.score = cur_score ;
02596 
02597     /* add column to hash table, for supercolumn detection */
02598     hash %= n_col + 1 ;
02599 
02600     DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
02601     ASSERT (((Int) hash) <= n_col) ;
02602 
02603     head_column = head [hash] ;
02604     if (head_column > EMPTY)
02605     {
02606         /* degree list "hash" is non-empty, use prev (shared3) of */
02607         /* first column in degree list as head of hash bucket */
02608         first_col = Col [head_column].shared3.headhash ;
02609         Col [head_column].shared3.headhash = col ;
02610     }
02611     else
02612     {
02613         /* degree list "hash" is empty, use head as hash bucket */
02614         first_col = - (head_column + 2) ;
02615         head [hash] = - (col + 2) ;
02616     }
02617     Col [col].shared4.hash_next = first_col ;
02618 
02619     /* save hash function in Col [col].shared3.hash */
02620     Col [col].shared3.hash = (Int) hash ;
02621     ASSERT (COL_IS_ALIVE (col)) ;
02622       }
02623   }
02624 
02625   /* The approximate external column degree is now computed.  */
02626 
02627   /* === Supercolumn detection ======================================== */
02628 
02629   DEBUG3 (("** Supercolumn detection phase. **\n")) ;
02630 
02631   detect_super_cols (
02632 
02633 #ifndef NDEBUG
02634     n_col, Row,
02635 #endif /* NDEBUG */
02636 
02637     Col, A, head, pivot_row_start, pivot_row_length) ;
02638 
02639   /* === Kill the pivotal column ====================================== */
02640 
02641   KILL_PRINCIPAL_COL (pivot_col) ;
02642 
02643   /* === Clear mark =================================================== */
02644 
02645   tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
02646 
02647 #ifndef NDEBUG
02648   DEBUG3 (("check3\n")) ;
02649   debug_mark (n_row, Row, tag_mark, max_mark) ;
02650 #endif /* NDEBUG */
02651 
02652   /* === Finalize the new pivot row, and column scores ================ */
02653 
02654   DEBUG3 (("** Finalize scores phase. **\n")) ;
02655 
02656   /* for each column in pivot row */
02657   rp = &A [pivot_row_start] ;
02658   /* compact the pivot row */
02659   new_rp = rp ;
02660   rp_end = rp + pivot_row_length ;
02661   while (rp < rp_end)
02662   {
02663       col = *rp++ ;
02664       /* skip dead columns */
02665       if (COL_IS_DEAD (col))
02666       {
02667     continue ;
02668       }
02669       *new_rp++ = col ;
02670       /* add new pivot row to column */
02671       A [Col [col].start + (Col [col].length++)] = pivot_row ;
02672 
02673       /* retrieve score so far and add on pivot row's degree. */
02674       /* (we wait until here for this in case the pivot */
02675       /* row's degree was reduced due to mass elimination). */
02676       cur_score = Col [col].shared2.score + pivot_row_degree ;
02677 
02678       /* calculate the max possible score as the number of */
02679       /* external columns minus the 'k' value minus the */
02680       /* columns thickness */
02681       max_score = n_col - k - Col [col].shared1.thickness ;
02682 
02683       /* make the score the external degree of the union-of-rows */
02684       cur_score -= Col [col].shared1.thickness ;
02685 
02686       /* make sure score is less or equal than the max score */
02687       cur_score = MIN (cur_score, max_score) ;
02688       ASSERT (cur_score >= 0) ;
02689 
02690       /* store updated score */
02691       Col [col].shared2.score = cur_score ;
02692 
02693       /* === Place column back in degree list ========================= */
02694 
02695       ASSERT (min_score >= 0) ;
02696       ASSERT (min_score <= n_col) ;
02697       ASSERT (cur_score >= 0) ;
02698       ASSERT (cur_score <= n_col) ;
02699       ASSERT (head [cur_score] >= EMPTY) ;
02700       next_col = head [cur_score] ;
02701       Col [col].shared4.degree_next = next_col ;
02702       Col [col].shared3.prev = EMPTY ;
02703       if (next_col != EMPTY)
02704       {
02705     Col [next_col].shared3.prev = col ;
02706       }
02707       head [cur_score] = col ;
02708 
02709       /* see if this score is less than current min */
02710       min_score = MIN (min_score, cur_score) ;
02711 
02712   }
02713 
02714 #ifndef NDEBUG
02715   debug_deg_lists (n_row, n_col, Row, Col, head,
02716     min_score, n_col2-k, max_deg) ;
02717 #endif /* NDEBUG */
02718 
02719   /* === Resurrect the new pivot row ================================== */
02720 
02721   if (pivot_row_degree > 0)
02722   {
02723       /* update pivot row length to reflect any cols that were killed */
02724       /* during super-col detection and mass elimination */
02725       Row [pivot_row].start  = pivot_row_start ;
02726       Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
02727       ASSERT (Row [pivot_row].length > 0) ;
02728       Row [pivot_row].shared1.degree = pivot_row_degree ;
02729       Row [pivot_row].shared2.mark = 0 ;
02730       /* pivot row is no longer dead */
02731 
02732       DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
02733       pivot_row, pivot_row_degree)) ;
02734   }
02735     }
02736 
02737     /* === All principal columns have now been ordered ====================== */
02738 
02739     return (ngarbage) ;
02740 }
02741 
02742 
02743 /* ========================================================================== */
02744 /* === order_children ======================================================= */
02745 /* ========================================================================== */
02746 
02747 /*
02748     The find_ordering routine has ordered all of the principal columns (the
02749     representatives of the supercolumns).  The non-principal columns have not
02750     yet been ordered.  This routine orders those columns by walking up the
02751     parent tree (a column is a child of the column which absorbed it).  The
02752     final permutation vector is then placed in p [0 ... n_col-1], with p [0]
02753     being the first column, and p [n_col-1] being the last.  It doesn't look
02754     like it at first glance, but be assured that this routine takes time linear
02755     in the number of columns.  Although not immediately obvious, the time
02756     taken by this routine is O (n_col), that is, linear in the number of
02757     columns.  Not user-callable.
02758 */
02759 
02760 PRIVATE void order_children
02761 (
02762     /* === Parameters ======================================================= */
02763 
02764     Int n_col,      /* number of columns of A */
02765     Colamd_Col Col [],    /* of size n_col+1 */
02766     Int p []      /* p [0 ... n_col-1] is the column permutation*/
02767 )
02768 {
02769     /* === Local variables ================================================== */
02770 
02771     Int i ;     /* loop counter for all columns */
02772     Int c ;     /* column index */
02773     Int parent ;    /* index of column's parent */
02774     Int order ;     /* column's order */
02775 
02776     /* === Order each non-principal column ================================== */
02777 
02778     for (i = 0 ; i < n_col ; i++)
02779     {
02780   /* find an un-ordered non-principal column */
02781   ASSERT (COL_IS_DEAD (i)) ;
02782   if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
02783   {
02784       parent = i ;
02785       /* once found, find its principal parent */
02786       do
02787       {
02788     parent = Col [parent].shared1.parent ;
02789       } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
02790 
02791       /* now, order all un-ordered non-principal columns along path */
02792       /* to this parent.  collapse tree at the same time */
02793       c = i ;
02794       /* get order of parent */
02795       order = Col [parent].shared2.order ;
02796 
02797       do
02798       {
02799     ASSERT (Col [c].shared2.order == EMPTY) ;
02800 
02801     /* order this column */
02802     Col [c].shared2.order = order++ ;
02803     /* collaps tree */
02804     Col [c].shared1.parent = parent ;
02805 
02806     /* get immediate parent of this column */
02807     c = Col [c].shared1.parent ;
02808 
02809     /* continue until we hit an ordered column.  There are */
02810     /* guarranteed not to be anymore unordered columns */
02811     /* above an ordered column */
02812       } while (Col [c].shared2.order == EMPTY) ;
02813 
02814       /* re-order the super_col parent to largest order for this group */
02815       Col [parent].shared2.order = order ;
02816   }
02817     }
02818 
02819     /* === Generate the permutation ========================================= */
02820 
02821     for (c = 0 ; c < n_col ; c++)
02822     {
02823   p [Col [c].shared2.order] = c ;
02824     }
02825 }
02826 
02827 
02828 /* ========================================================================== */
02829 /* === detect_super_cols ==================================================== */
02830 /* ========================================================================== */
02831 
02832 /*
02833     Detects supercolumns by finding matches between columns in the hash buckets.
02834     Check amongst columns in the set A [row_start ... row_start + row_length-1].
02835     The columns under consideration are currently *not* in the degree lists,
02836     and have already been placed in the hash buckets.
02837 
02838     The hash bucket for columns whose hash function is equal to h is stored
02839     as follows:
02840 
02841   if head [h] is >= 0, then head [h] contains a degree list, so:
02842 
02843     head [h] is the first column in degree bucket h.
02844     Col [head [h]].headhash gives the first column in hash bucket h.
02845 
02846   otherwise, the degree list is empty, and:
02847 
02848     -(head [h] + 2) is the first column in hash bucket h.
02849 
02850     For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
02851     column" pointer.  Col [c].shared3.hash is used instead as the hash number
02852     for that column.  The value of Col [c].shared4.hash_next is the next column
02853     in the same hash bucket.
02854 
02855     Assuming no, or "few" hash collisions, the time taken by this routine is
02856     linear in the sum of the sizes (lengths) of each column whose score has
02857     just been computed in the approximate degree computation.
02858     Not user-callable.
02859 */
02860 
02861 PRIVATE void detect_super_cols
02862 (
02863     /* === Parameters ======================================================= */
02864 
02865 #ifndef NDEBUG
02866     /* these two parameters are only needed when debugging is enabled: */
02867     Int n_col,      /* number of columns of A */
02868     Colamd_Row Row [],    /* of size n_row+1 */
02869 #endif /* NDEBUG */
02870 
02871     Colamd_Col Col [],    /* of size n_col+1 */
02872     Int A [],     /* row indices of A */
02873     Int head [],    /* head of degree lists and hash buckets */
02874     Int row_start,    /* pointer to set of columns to check */
02875     Int row_length    /* number of columns to check */
02876 )
02877 {
02878     /* === Local variables ================================================== */
02879 
02880     Int hash ;      /* hash value for a column */
02881     Int *rp ;     /* pointer to a row */
02882     Int c ;     /* a column index */
02883     Int super_c ;   /* column index of the column to absorb into */
02884     Int *cp1 ;      /* column pointer for column super_c */
02885     Int *cp2 ;      /* column pointer for column c */
02886     Int length ;    /* length of column super_c */
02887     Int prev_c ;    /* column preceding c in hash bucket */
02888     Int i ;     /* loop counter */
02889     Int *rp_end ;   /* pointer to the end of the row */
02890     Int col ;     /* a column index in the row to check */
02891     Int head_column ;   /* first column in hash bucket or degree list */
02892     Int first_col ;   /* first column in hash bucket */
02893 
02894     /* === Consider each column in the row ================================== */
02895 
02896     rp = &A [row_start] ;
02897     rp_end = rp + row_length ;
02898     while (rp < rp_end)
02899     {
02900   col = *rp++ ;
02901   if (COL_IS_DEAD (col))
02902   {
02903       continue ;
02904   }
02905 
02906   /* get hash number for this column */
02907   hash = Col [col].shared3.hash ;
02908   ASSERT (hash <= n_col) ;
02909 
02910   /* === Get the first column in this hash bucket ===================== */
02911 
02912   head_column = head [hash] ;
02913   if (head_column > EMPTY)
02914   {
02915       first_col = Col [head_column].shared3.headhash ;
02916   }
02917   else
02918   {
02919       first_col = - (head_column + 2) ;
02920   }
02921 
02922   /* === Consider each column in the hash bucket ====================== */
02923 
02924   for (super_c = first_col ; super_c != EMPTY ;
02925       super_c = Col [super_c].shared4.hash_next)
02926   {
02927       ASSERT (COL_IS_ALIVE (super_c)) ;
02928       ASSERT (Col [super_c].shared3.hash == hash) ;
02929       length = Col [super_c].length ;
02930 
02931       /* prev_c is the column preceding column c in the hash bucket */
02932       prev_c = super_c ;
02933 
02934       /* === Compare super_c with all columns after it ================ */
02935 
02936       for (c = Col [super_c].shared4.hash_next ;
02937      c != EMPTY ; c = Col [c].shared4.hash_next)
02938       {
02939     ASSERT (c != super_c) ;
02940     ASSERT (COL_IS_ALIVE (c)) ;
02941     ASSERT (Col [c].shared3.hash == hash) ;
02942 
02943     /* not identical if lengths or scores are different */
02944     if (Col [c].length != length ||
02945         Col [c].shared2.score != Col [super_c].shared2.score)
02946     {
02947         prev_c = c ;
02948         continue ;
02949     }
02950 
02951     /* compare the two columns */
02952     cp1 = &A [Col [super_c].start] ;
02953     cp2 = &A [Col [c].start] ;
02954 
02955     for (i = 0 ; i < length ; i++)
02956     {
02957         /* the columns are "clean" (no dead rows) */
02958         ASSERT (ROW_IS_ALIVE (*cp1))  ;
02959         ASSERT (ROW_IS_ALIVE (*cp2))  ;
02960         /* row indices will same order for both supercols, */
02961         /* no gather scatter nessasary */
02962         if (*cp1++ != *cp2++)
02963         {
02964       break ;
02965         }
02966     }
02967 
02968     /* the two columns are different if the for-loop "broke" */
02969     if (i != length)
02970     {
02971         prev_c = c ;
02972         continue ;
02973     }
02974 
02975     /* === Got it!  two columns are identical =================== */
02976 
02977     ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
02978 
02979     Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
02980     Col [c].shared1.parent = super_c ;
02981     KILL_NON_PRINCIPAL_COL (c) ;
02982     /* order c later, in order_children() */
02983     Col [c].shared2.order = EMPTY ;
02984     /* remove c from hash bucket */
02985     Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
02986       }
02987   }
02988 
02989   /* === Empty this hash bucket ======================================= */
02990 
02991   if (head_column > EMPTY)
02992   {
02993       /* corresponding degree list "hash" is not empty */
02994       Col [head_column].shared3.headhash = EMPTY ;
02995   }
02996   else
02997   {
02998       /* corresponding degree list "hash" is empty */
02999       head [hash] = EMPTY ;
03000   }
03001     }
03002 }
03003 
03004 
03005 /* ========================================================================== */
03006 /* === garbage_collection =================================================== */
03007 /* ========================================================================== */
03008 
03009 /*
03010     Defragments and compacts columns and rows in the workspace A.  Used when
03011     all avaliable memory has been used while performing row merging.  Returns
03012     the index of the first free position in A, after garbage collection.  The
03013     time taken by this routine is linear is the size of the array A, which is
03014     itself linear in the number of nonzeros in the input matrix.
03015     Not user-callable.
03016 */
03017 
03018 PRIVATE Int garbage_collection  /* returns the new value of pfree */
03019 (
03020     /* === Parameters ======================================================= */
03021 
03022     Int n_row,      /* number of rows */
03023     Int n_col,      /* number of columns */
03024     Colamd_Row Row [],    /* row info */
03025     Colamd_Col Col [],    /* column info */
03026     Int A [],     /* A [0 ... Alen-1] holds the matrix */
03027     Int *pfree      /* &A [0] ... pfree is in use */
03028 )
03029 {
03030     /* === Local variables ================================================== */
03031 
03032     Int *psrc ;     /* source pointer */
03033     Int *pdest ;    /* destination pointer */
03034     Int j ;     /* counter */
03035     Int r ;     /* a row index */
03036     Int c ;     /* a column index */
03037     Int length ;    /* length of a row or column */
03038 
03039 #ifndef NDEBUG
03040     Int debug_rows ;
03041     DEBUG2 (("Defrag..\n")) ;
03042     for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
03043     debug_rows = 0 ;
03044 #endif /* NDEBUG */
03045 
03046     /* === Defragment the columns =========================================== */
03047 
03048     pdest = &A[0] ;
03049     for (c = 0 ; c < n_col ; c++)
03050     {
03051   if (COL_IS_ALIVE (c))
03052   {
03053       psrc = &A [Col [c].start] ;
03054 
03055       /* move and compact the column */
03056       ASSERT (pdest <= psrc) ;
03057       Col [c].start = (Int) (pdest - &A [0]) ;
03058       length = Col [c].length ;
03059       for (j = 0 ; j < length ; j++)
03060       {
03061     r = *psrc++ ;
03062     if (ROW_IS_ALIVE (r))
03063     {
03064         *pdest++ = r ;
03065     }
03066       }
03067       Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
03068   }
03069     }
03070 
03071     /* === Prepare to defragment the rows =================================== */
03072 
03073     for (r = 0 ; r < n_row ; r++)
03074     {
03075   if (ROW_IS_DEAD (r) || (Row [r].length == 0))
03076   {
03077       /* This row is already dead, or is of zero length.  Cannot compact
03078        * a row of zero length, so kill it.  NOTE: in the current version,
03079        * there are no zero-length live rows.  Kill the row (for the first
03080        * time, or again) just to be safe. */
03081       KILL_ROW (r) ;
03082   }
03083   else
03084   {
03085       /* save first column index in Row [r].shared2.first_column */
03086       psrc = &A [Row [r].start] ;
03087       Row [r].shared2.first_column = *psrc ;
03088       ASSERT (ROW_IS_ALIVE (r)) ;
03089       /* flag the start of the row with the one's complement of row */
03090       *psrc = ONES_COMPLEMENT (r) ;
03091 #ifndef NDEBUG
03092       debug_rows++ ;
03093 #endif /* NDEBUG */
03094   }
03095     }
03096 
03097     /* === Defragment the rows ============================================== */
03098 
03099     psrc = pdest ;
03100     while (psrc < pfree)
03101     {
03102   /* find a negative number ... the start of a row */
03103   if (*psrc++ < 0)
03104   {
03105       psrc-- ;
03106       /* get the row index */
03107       r = ONES_COMPLEMENT (*psrc) ;
03108       ASSERT (r >= 0 && r < n_row) ;
03109       /* restore first column index */
03110       *psrc = Row [r].shared2.first_column ;
03111       ASSERT (ROW_IS_ALIVE (r)) ;
03112       ASSERT (Row [r].length > 0) ;
03113       /* move and compact the row */
03114       ASSERT (pdest <= psrc) ;
03115       Row [r].start = (Int) (pdest - &A [0]) ;
03116       length = Row [r].length ;
03117       for (j = 0 ; j < length ; j++)
03118       {
03119     c = *psrc++ ;
03120     if (COL_IS_ALIVE (c))
03121     {
03122         *pdest++ = c ;
03123     }
03124       }
03125       Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
03126       ASSERT (Row [r].length > 0) ;
03127 #ifndef NDEBUG
03128       debug_rows-- ;
03129 #endif /* NDEBUG */
03130   }
03131     }
03132     /* ensure we found all the rows */
03133     ASSERT (debug_rows == 0) ;
03134 
03135     /* === Return the new value of pfree ==================================== */
03136 
03137     return ((Int) (pdest - &A [0])) ;
03138 }
03139 
03140 
03141 /* ========================================================================== */
03142 /* === clear_mark =========================================================== */
03143 /* ========================================================================== */
03144 
03145 /*
03146     Clears the Row [].shared2.mark array, and returns the new tag_mark.
03147     Return value is the new tag_mark.  Not user-callable.
03148 */
03149 
03150 PRIVATE Int clear_mark  /* return the new value for tag_mark */
03151 (
03152     /* === Parameters ======================================================= */
03153 
03154     Int tag_mark, /* new value of tag_mark */
03155     Int max_mark, /* max allowed value of tag_mark */
03156 
03157     Int n_row,    /* number of rows in A */
03158     Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
03159 )
03160 {
03161     /* === Local variables ================================================== */
03162 
03163     Int r ;
03164 
03165     if (tag_mark <= 0 || tag_mark >= max_mark)
03166     {
03167   for (r = 0 ; r < n_row ; r++)
03168   {
03169       if (ROW_IS_ALIVE (r))
03170       {
03171     Row [r].shared2.mark = 0 ;
03172       }
03173   }
03174   tag_mark = 1 ;
03175     }
03176 
03177     return (tag_mark) ;
03178 }
03179 
03180 
03181 /* ========================================================================== */
03182 /* === print_report ========================================================= */
03183 /* ========================================================================== */
03184 
03185 PRIVATE void print_report
03186 (
03187     char *method,
03188     Int stats [COLAMD_STATS]
03189 )
03190 {
03191 
03192     Int i1, i2, i3 ;
03193 
03194     PRINTF (("\n%s version %d.%d, %s: ", method,
03195       COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
03196 
03197     if (!stats)
03198     {
03199       PRINTF (("No statistics available.\n")) ;
03200   return ;
03201     }
03202 
03203     i1 = stats [COLAMD_INFO1] ;
03204     i2 = stats [COLAMD_INFO2] ;
03205     i3 = stats [COLAMD_INFO3] ;
03206 
03207     if (stats [COLAMD_STATUS] >= 0)
03208     {
03209       PRINTF (("OK.  ")) ;
03210     }
03211     else
03212     {
03213       PRINTF (("ERROR.  ")) ;
03214     }
03215 
03216     switch (stats [COLAMD_STATUS])
03217     {
03218 
03219   case COLAMD_OK_BUT_JUMBLED:
03220 
03221       PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
03222 
03223       PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
03224       method, i3)) ;
03225 
03226       PRINTF(("%s: last seen duplicate or out-of-order row index:   %d\n",
03227       method, INDEX (i2))) ;
03228 
03229       PRINTF(("%s: last seen in column:                             %d",
03230       method, INDEX (i1))) ;
03231 
03232       /* no break - fall through to next case instead */
03233 
03234   case COLAMD_OK:
03235 
03236       PRINTF(("\n")) ;
03237 
03238       PRINTF(("%s: number of dense or empty rows ignored:           %d\n",
03239       method, stats [COLAMD_DENSE_ROW])) ;
03240 
03241       PRINTF(("%s: number of dense or empty columns ignored:        %d\n",
03242       method, stats [COLAMD_DENSE_COL])) ;
03243 
03244       PRINTF(("%s: number of garbage collections performed:         %d\n",
03245       method, stats [COLAMD_DEFRAG_COUNT])) ;
03246       break ;
03247 
03248   case COLAMD_ERROR_A_not_present:
03249 
03250       PRINTF(("Array A (row indices of matrix) not present.\n")) ;
03251       break ;
03252 
03253   case COLAMD_ERROR_p_not_present:
03254 
03255       PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
03256       break ;
03257 
03258   case COLAMD_ERROR_nrow_negative:
03259 
03260       PRINTF(("Invalid number of rows (%d).\n", i1)) ;
03261       break ;
03262 
03263   case COLAMD_ERROR_ncol_negative:
03264 
03265       PRINTF(("Invalid number of columns (%d).\n", i1)) ;
03266       break ;
03267 
03268   case COLAMD_ERROR_nnz_negative:
03269 
03270       PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
03271       break ;
03272 
03273   case COLAMD_ERROR_p0_nonzero:
03274 
03275       PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
03276       break ;
03277 
03278   case COLAMD_ERROR_A_too_small:
03279 
03280       PRINTF(("Array A too small.\n")) ;
03281       PRINTF(("        Need Alen >= %d, but given only Alen = %d.\n",
03282       i1, i2)) ;
03283       break ;
03284 
03285   case COLAMD_ERROR_col_length_negative:
03286 
03287       PRINTF
03288       (("Column %d has a negative number of nonzero entries (%d).\n",
03289       INDEX (i1), i2)) ;
03290       break ;
03291 
03292   case COLAMD_ERROR_row_index_out_of_bounds:
03293 
03294       PRINTF
03295       (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
03296       INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
03297       break ;
03298 
03299   case COLAMD_ERROR_out_of_memory:
03300 
03301       PRINTF(("Out of memory.\n")) ;
03302       break ;
03303 
03304   /* v2.4: internal-error case deleted */
03305     }
03306 }
03307 
03308 
03309 
03310 
03311 /* ========================================================================== */
03312 /* === colamd debugging routines ============================================ */
03313 /* ========================================================================== */
03314 
03315 /* When debugging is disabled, the remainder of this file is ignored. */
03316 
03317 #ifndef NDEBUG
03318 
03319 
03320 /* ========================================================================== */
03321 /* === debug_structures ===================================================== */
03322 /* ========================================================================== */
03323 
03324 /*
03325     At this point, all empty rows and columns are dead.  All live columns
03326     are "clean" (containing no dead rows) and simplicial (no supercolumns
03327     yet).  Rows may contain dead columns, but all live rows contain at
03328     least one live column.
03329 */
03330 
03331 PRIVATE void debug_structures
03332 (
03333     /* === Parameters ======================================================= */
03334 
03335     Int n_row,
03336     Int n_col,
03337     Colamd_Row Row [],
03338     Colamd_Col Col [],
03339     Int A [],
03340     Int n_col2
03341 )
03342 {
03343     /* === Local variables ================================================== */
03344 
03345     Int i ;
03346     Int c ;
03347     Int *cp ;
03348     Int *cp_end ;
03349     Int len ;
03350     Int score ;
03351     Int r ;
03352     Int *rp ;
03353     Int *rp_end ;
03354     Int deg ;
03355 
03356     /* === Check A, Row, and Col ============================================ */
03357 
03358     for (c = 0 ; c < n_col ; c++)
03359     {
03360   if (COL_IS_ALIVE (c))
03361   {
03362       len = Col [c].length ;
03363       score = Col [c].shared2.score ;
03364       DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
03365       ASSERT (len > 0) ;
03366       ASSERT (score >= 0) ;
03367       ASSERT (Col [c].shared1.thickness == 1) ;
03368       cp = &A [Col [c].start] ;
03369       cp_end = cp + len ;
03370       while (cp < cp_end)
03371       {
03372     r = *cp++ ;
03373     ASSERT (ROW_IS_ALIVE (r)) ;
03374       }
03375   }
03376   else
03377   {
03378       i = Col [c].shared2.order ;
03379       ASSERT (i >= n_col2 && i < n_col) ;
03380   }
03381     }
03382 
03383     for (r = 0 ; r < n_row ; r++)
03384     {
03385   if (ROW_IS_ALIVE (r))
03386   {
03387       i = 0 ;
03388       len = Row [r].length ;
03389       deg = Row [r].shared1.degree ;
03390       ASSERT (len > 0) ;
03391       ASSERT (deg > 0) ;
03392       rp = &A [Row [r].start] ;
03393       rp_end = rp + len ;
03394       while (rp < rp_end)
03395       {
03396     c = *rp++ ;
03397     if (COL_IS_ALIVE (c))
03398     {
03399         i++ ;
03400     }
03401       }
03402       ASSERT (i > 0) ;
03403   }
03404     }
03405 }
03406 
03407 
03408 /* ========================================================================== */
03409 /* === debug_deg_lists ====================================================== */
03410 /* ========================================================================== */
03411 
03412 /*
03413     Prints the contents of the degree lists.  Counts the number of columns
03414     in the degree list and compares it to the total it should have.  Also
03415     checks the row degrees.
03416 */
03417 
03418 PRIVATE void debug_deg_lists
03419 (
03420     /* === Parameters ======================================================= */
03421 
03422     Int n_row,
03423     Int n_col,
03424     Colamd_Row Row [],
03425     Colamd_Col Col [],
03426     Int head [],
03427     Int min_score,
03428     Int should,
03429     Int max_deg
03430 )
03431 {
03432     /* === Local variables ================================================== */
03433 
03434     Int deg ;
03435     Int col ;
03436     Int have ;
03437     Int row ;
03438 
03439     /* === Check the degree lists =========================================== */
03440 
03441     if (n_col > 10000 && colamd_debug <= 0)
03442     {
03443   return ;
03444     }
03445     have = 0 ;
03446     DEBUG4 (("Degree lists: %d\n", min_score)) ;
03447     for (deg = 0 ; deg <= n_col ; deg++)
03448     {
03449   col = head [deg] ;
03450   if (col == EMPTY)
03451   {
03452       continue ;
03453   }
03454   DEBUG4 (("%d:", deg)) ;
03455   while (col != EMPTY)
03456   {
03457       DEBUG4 ((" %d", col)) ;
03458       have += Col [col].shared1.thickness ;
03459       ASSERT (COL_IS_ALIVE (col)) ;
03460       col = Col [col].shared4.degree_next ;
03461   }
03462   DEBUG4 (("\n")) ;
03463     }
03464     DEBUG4 (("should %d have %d\n", should, have)) ;
03465     ASSERT (should == have) ;
03466 
03467     /* === Check the row degrees ============================================ */
03468 
03469     if (n_row > 10000 && colamd_debug <= 0)
03470     {
03471   return ;
03472     }
03473     for (row = 0 ; row < n_row ; row++)
03474     {
03475   if (ROW_IS_ALIVE (row))
03476   {
03477       ASSERT (Row [row].shared1.degree <= max_deg) ;
03478   }
03479     }
03480 }
03481 
03482 
03483 /* ========================================================================== */
03484 /* === debug_mark =========================================================== */
03485 /* ========================================================================== */
03486 
03487 /*
03488     Ensures that the tag_mark is less that the maximum and also ensures that
03489     each entry in the mark array is less than the tag mark.
03490 */
03491 
03492 PRIVATE void debug_mark
03493 (
03494     /* === Parameters ======================================================= */
03495 
03496     Int n_row,
03497     Colamd_Row Row [],
03498     Int tag_mark,
03499     Int max_mark
03500 )
03501 {
03502     /* === Local variables ================================================== */
03503 
03504     Int r ;
03505 
03506     /* === Check the Row marks ============================================== */
03507 
03508     ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
03509     if (n_row > 10000 && colamd_debug <= 0)
03510     {
03511   return ;
03512     }
03513     for (r = 0 ; r < n_row ; r++)
03514     {
03515   ASSERT (Row [r].shared2.mark < tag_mark) ;
03516     }
03517 }
03518 
03519 
03520 /* ========================================================================== */
03521 /* === debug_matrix ========================================================= */
03522 /* ========================================================================== */
03523 
03524 /*
03525     Prints out the contents of the columns and the rows.
03526 */
03527 
03528 PRIVATE void debug_matrix
03529 (
03530     /* === Parameters ======================================================= */
03531 
03532     Int n_row,
03533     Int n_col,
03534     Colamd_Row Row [],
03535     Colamd_Col Col [],
03536     Int A []
03537 )
03538 {
03539     /* === Local variables ================================================== */
03540 
03541     Int r ;
03542     Int c ;
03543     Int *rp ;
03544     Int *rp_end ;
03545     Int *cp ;
03546     Int *cp_end ;
03547 
03548     /* === Dump the rows and columns of the matrix ========================== */
03549 
03550     if (colamd_debug < 3)
03551     {
03552   return ;
03553     }
03554     DEBUG3 (("DUMP MATRIX:\n")) ;
03555     for (r = 0 ; r < n_row ; r++)
03556     {
03557   DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
03558   if (ROW_IS_DEAD (r))
03559   {
03560       continue ;
03561   }
03562   DEBUG3 (("start %d length %d degree %d\n",
03563     Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
03564   rp = &A [Row [r].start] ;
03565   rp_end = rp + Row [r].length ;
03566   while (rp < rp_end)
03567   {
03568       c = *rp++ ;
03569       DEBUG4 (("  %d col %d\n", COL_IS_ALIVE (c), c)) ;
03570   }
03571     }
03572 
03573     for (c = 0 ; c < n_col ; c++)
03574     {
03575   DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
03576   if (COL_IS_DEAD (c))
03577   {
03578       continue ;
03579   }
03580   DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
03581     Col [c].start, Col [c].length,
03582     Col [c].shared1.thickness, Col [c].shared2.score)) ;
03583   cp = &A [Col [c].start] ;
03584   cp_end = cp + Col [c].length ;
03585   while (cp < cp_end)
03586   {
03587       r = *cp++ ;
03588       DEBUG4 (("  %d row %d\n", ROW_IS_ALIVE (r), r)) ;
03589   }
03590     }
03591 }
03592 
03593 PRIVATE void colamd_get_debug
03594 (
03595     char *method
03596 )
03597 {
03598     FILE *f ;
03599     colamd_debug = 0 ;    /* no debug printing */
03600     f = fopen ("debug", "r") ;
03601     if (f == (FILE *) NULL)
03602     {
03603   colamd_debug = 0 ;
03604     }
03605     else
03606     {
03607   fscanf (f, "%d", &colamd_debug) ;
03608   fclose (f) ;
03609     }
03610     DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
03611       method, colamd_debug)) ;
03612 }
03613 
03614 #endif /* NDEBUG */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines