Amesos Package Browser (Single Doxygen Collection) Development
amesos_colamd.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 #include "amesos_colamd.h"
00667 #include <limits.h>
00668 #include <math.h>
00669 
00670 #ifdef MATLAB_MEX_FILE
00671 #include "mex.h"
00672 #include "matrix.h"
00673 #endif /* MATLAB_MEX_FILE */
00674 
00675 #if !defined (NPRINT) || !defined (NDEBUG)
00676 #include <stdio.h>
00677 #endif
00678 
00679 #ifndef NULL
00680 #define NULL ((void *) 0)
00681 #endif
00682 
00683 /* ========================================================================== */
00684 /* === int or UF_long ======================================================= */
00685 /* ========================================================================== */
00686 
00687 /* define UF_long */
00688 #include "amesos_UFconfig.h"
00689 
00690 #ifdef DLONG
00691 
00692 #define Int UF_long
00693 #define ID  UF_long_id
00694 #define Int_MAX UF_long_max
00695 
00696 #define COLAMD_recommended amesos_colamd_l_recommended
00697 #define COLAMD_set_defaults amesos_colamd_l_set_defaults
00698 #define COLAMD_MAIN amesos_colamd_l
00699 #define SYMAMD_MAIN amesos_symamd_l
00700 #define COLAMD_report amesos_colamd_l_report
00701 #define SYMAMD_report amesos_symamd_l_report
00702 
00703 #else
00704 
00705 #define Int int
00706 #define ID "%d"
00707 #define Int_MAX INT_MAX
00708 
00709 #define COLAMD_recommended amesos_colamd_recommended
00710 #define COLAMD_set_defaults amesos_colamd_set_defaults
00711 #define COLAMD_MAIN amesos_colamd
00712 #define SYMAMD_MAIN amesos_symamd
00713 #define COLAMD_report amesos_colamd_report
00714 #define SYMAMD_report amesos_symamd_report
00715 
00716 #endif
00717 
00718 /* ========================================================================== */
00719 /* === Row and Column structures ============================================ */
00720 /* ========================================================================== */
00721 
00722 /* User code that makes use of the colamd/symamd routines need not directly */
00723 /* reference these structures.  They are used only for colamd_recommended. */
00724 
00725 typedef struct Colamd_Col_struct
00726 {
00727     Int start ;   /* index for A of first row in this column, or DEAD */
00728       /* if column is dead */
00729     Int length ;  /* number of rows in this column */
00730     union
00731     {
00732   Int thickness ; /* number of original columns represented by this */
00733       /* col, if the column is alive */
00734   Int parent ;  /* parent in parent tree super-column structure, if */
00735       /* the column is dead */
00736     } shared1 ;
00737     union
00738     {
00739   Int score ; /* the score used to maintain heap, if col is alive */
00740   Int order ; /* pivot ordering of this column, if col is dead */
00741     } shared2 ;
00742     union
00743     {
00744   Int headhash ;  /* head of a hash bucket, if col is at the head of */
00745       /* a degree list */
00746   Int hash ;  /* hash value, if col is not in a degree list */
00747   Int prev ;  /* previous column in degree list, if col is in a */
00748       /* degree list (but not at the head of a degree list) */
00749     } shared3 ;
00750     union
00751     {
00752   Int degree_next ; /* next column, if col is in a degree list */
00753   Int hash_next ;   /* next column, if col is in a hash list */
00754     } shared4 ;
00755 
00756 } Colamd_Col ;
00757 
00758 typedef struct Colamd_Row_struct
00759 {
00760     Int start ;   /* index for A of first col in this row */
00761     Int length ;  /* number of principal columns in this row */
00762     union
00763     {
00764   Int degree ;  /* number of principal & non-principal columns in row */
00765   Int p ;   /* used as a row pointer in init_rows_cols () */
00766     } shared1 ;
00767     union
00768     {
00769   Int mark ;  /* for computing set differences and marking dead rows*/
00770   Int first_column ;/* first column in row (used in garbage collection) */
00771     } shared2 ;
00772 
00773 } Colamd_Row ;
00774 
00775 /* ========================================================================== */
00776 /* === Definitions ========================================================== */
00777 /* ========================================================================== */
00778 
00779 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
00780 #define PUBLIC
00781 #define PRIVATE static
00782 
00783 #define DENSE_DEGREE(alpha,n) \
00784     ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
00785 
00786 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
00787 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
00788 
00789 #define ONES_COMPLEMENT(r) (-(r)-1)
00790 
00791 /* -------------------------------------------------------------------------- */
00792 /* Change for version 2.1:  define TRUE and FALSE only if not yet defined */  
00793 /* -------------------------------------------------------------------------- */
00794 
00795 #ifndef TRUE
00796 #define TRUE (1)
00797 #endif
00798 
00799 #ifndef FALSE
00800 #define FALSE (0)
00801 #endif
00802 
00803 /* -------------------------------------------------------------------------- */
00804 
00805 #define EMPTY (-1)
00806 
00807 /* Row and column status */
00808 #define ALIVE (0)
00809 #define DEAD  (-1)
00810 
00811 /* Column status */
00812 #define DEAD_PRINCIPAL    (-1)
00813 #define DEAD_NON_PRINCIPAL  (-2)
00814 
00815 /* Macros for row and column status update and checking. */
00816 #define ROW_IS_DEAD(r)      ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
00817 #define ROW_IS_MARKED_DEAD(row_mark)  (row_mark < ALIVE)
00818 #define ROW_IS_ALIVE(r)     (Row [r].shared2.mark >= ALIVE)
00819 #define COL_IS_DEAD(c)      (Col [c].start < ALIVE)
00820 #define COL_IS_ALIVE(c)     (Col [c].start >= ALIVE)
00821 #define COL_IS_DEAD_PRINCIPAL(c)  (Col [c].start == DEAD_PRINCIPAL)
00822 #define KILL_ROW(r)     { Row [r].shared2.mark = DEAD ; }
00823 #define KILL_PRINCIPAL_COL(c)   { Col [c].start = DEAD_PRINCIPAL ; }
00824 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
00825 
00826 /* ========================================================================== */
00827 /* === Colamd reporting mechanism =========================================== */
00828 /* ========================================================================== */
00829 
00830 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
00831 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
00832 #define INDEX(i) ((i)+1)
00833 #else
00834 /* In C, matrices are 0-based and indices are reported as such in *_report */
00835 #define INDEX(i) (i)
00836 #endif
00837 
00838 /* All output goes through the PRINTF macro.  */
00839 #define PRINTF(params) { if (amesos_colamd_printf != NULL) (void) amesos_colamd_printf params ; }
00840 
00841 /* ========================================================================== */
00842 /* === Prototypes of PRIVATE routines ======================================= */
00843 /* ========================================================================== */
00844 
00845 PRIVATE Int init_rows_cols
00846 (
00847     Int n_row,
00848     Int n_col,
00849     Colamd_Row Row [],
00850     Colamd_Col Col [],
00851     Int A [],
00852     Int p [],
00853     Int stats [COLAMD_STATS]
00854 ) ;
00855 
00856 PRIVATE void init_scoring
00857 (
00858     Int n_row,
00859     Int n_col,
00860     Colamd_Row Row [],
00861     Colamd_Col Col [],
00862     Int A [],
00863     Int head [],
00864     double knobs [COLAMD_KNOBS],
00865     Int *p_n_row2,
00866     Int *p_n_col2,
00867     Int *p_max_deg
00868 ) ;
00869 
00870 PRIVATE Int find_ordering
00871 (
00872     Int n_row,
00873     Int n_col,
00874     Int Alen,
00875     Colamd_Row Row [],
00876     Colamd_Col Col [],
00877     Int A [],
00878     Int head [],
00879     Int n_col2,
00880     Int max_deg,
00881     Int pfree,
00882     Int aggressive
00883 ) ;
00884 
00885 PRIVATE void order_children
00886 (
00887     Int n_col,
00888     Colamd_Col Col [],
00889     Int p []
00890 ) ;
00891 
00892 PRIVATE void detect_super_cols
00893 (
00894 
00895 #ifndef NDEBUG
00896     Int n_col,
00897     Colamd_Row Row [],
00898 #endif /* NDEBUG */
00899 
00900     Colamd_Col Col [],
00901     Int A [],
00902     Int head [],
00903     Int row_start,
00904     Int row_length
00905 ) ;
00906 
00907 PRIVATE Int garbage_collection
00908 (
00909     Int n_row,
00910     Int n_col,
00911     Colamd_Row Row [],
00912     Colamd_Col Col [],
00913     Int A [],
00914     Int *pfree
00915 ) ;
00916 
00917 PRIVATE Int clear_mark
00918 (
00919     Int tag_mark,
00920     Int max_mark,
00921     Int n_row,
00922     Colamd_Row Row []
00923 ) ;
00924 
00925 PRIVATE void print_report
00926 (
00927     char *method,
00928     Int stats [COLAMD_STATS]
00929 ) ;
00930 
00931 /* ========================================================================== */
00932 /* === Debugging prototypes and definitions ================================= */
00933 /* ========================================================================== */
00934 
00935 #ifndef NDEBUG
00936 
00937 #include <assert.h>
00938 
00939 /* colamd_debug is the *ONLY* global variable, and is only */
00940 /* present when debugging */
00941 
00942 PRIVATE Int colamd_debug = 0 ;  /* debug print level */
00943 
00944 #define DEBUG0(params) { PRINTF (params) ; }
00945 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
00946 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
00947 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
00948 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
00949 
00950 #ifdef MATLAB_MEX_FILE
00951 #define ASSERT(expression) (mxAssert ((expression), ""))
00952 #else
00953 #define ASSERT(expression) (assert (expression))
00954 #endif /* MATLAB_MEX_FILE */
00955 
00956 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */
00957 (
00958     char *method
00959 ) ;
00960 
00961 PRIVATE void debug_deg_lists
00962 (
00963     Int n_row,
00964     Int n_col,
00965     Colamd_Row Row [],
00966     Colamd_Col Col [],
00967     Int head [],
00968     Int min_score,
00969     Int should,
00970     Int max_deg
00971 ) ;
00972 
00973 PRIVATE void debug_mark
00974 (
00975     Int n_row,
00976     Colamd_Row Row [],
00977     Int tag_mark,
00978     Int max_mark
00979 ) ;
00980 
00981 PRIVATE void debug_matrix
00982 (
00983     Int n_row,
00984     Int n_col,
00985     Colamd_Row Row [],
00986     Colamd_Col Col [],
00987     Int A []
00988 ) ;
00989 
00990 PRIVATE void debug_structures
00991 (
00992     Int n_row,
00993     Int n_col,
00994     Colamd_Row Row [],
00995     Colamd_Col Col [],
00996     Int A [],
00997     Int n_col2
00998 ) ;
00999 
01000 #else /* NDEBUG */
01001 
01002 /* === No debugging ========================================================= */
01003 
01004 #define DEBUG0(params) ;
01005 #define DEBUG1(params) ;
01006 #define DEBUG2(params) ;
01007 #define DEBUG3(params) ;
01008 #define DEBUG4(params) ;
01009 
01010 #define ASSERT(expression)
01011 
01012 #endif /* NDEBUG */
01013 
01014 /* ========================================================================== */
01015 /* === USER-CALLABLE ROUTINES: ============================================== */
01016 /* ========================================================================== */
01017 
01018 /* ========================================================================== */
01019 /* === colamd_recommended =================================================== */
01020 /* ========================================================================== */
01021 
01022 /*
01023     The colamd_recommended routine returns the suggested size for Alen.  This
01024     value has been determined to provide good balance between the number of
01025     garbage collections and the memory requirements for colamd.  If any
01026     argument is negative, or if integer overflow occurs, a 0 is returned as an
01027     error condition.  2*nnz space is required for the row and column
01028     indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
01029     required for the Col and Row arrays, respectively, which are internal to
01030     colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the
01031     minimal amount of "elbow room", and nnz/5 more space is recommended for
01032     run time efficiency.
01033 
01034     Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
01035 
01036     This function is not needed when using symamd.
01037 */
01038 
01039 /* add two values of type size_t, and check for integer overflow */
01040 static size_t t_add (size_t a, size_t b, int *ok)
01041 {
01042     (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
01043     return ((*ok) ? (a + b) : 0) ;
01044 }
01045 
01046 /* compute a*k where k is a small integer, and check for integer overflow */
01047 static size_t t_mult (size_t a, size_t k, int *ok)
01048 {
01049     size_t i, s = 0 ;
01050     for (i = 0 ; i < k ; i++)
01051     {
01052   s = t_add (s, a, ok) ;
01053     }
01054     return (s) ;
01055 }
01056 
01057 /* size of the Col and Row structures */
01058 #define COLAMD_C(n_col,ok) \
01059     ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
01060 
01061 #define COLAMD_R(n_row,ok) \
01062     ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
01063 
01064 
01065 PUBLIC size_t COLAMD_recommended  /* returns recommended value of Alen. */
01066 (
01067     /* === Parameters ======================================================= */
01068 
01069     Int nnz,      /* number of nonzeros in A */
01070     Int n_row,      /* number of rows in A */
01071     Int n_col     /* number of columns in A */
01072 )
01073 {
01074     size_t s, c, r ;
01075     int ok = TRUE ;
01076     if (nnz < 0 || n_row < 0 || n_col < 0)
01077     {
01078   return (0) ;
01079     }
01080     s = t_mult (nnz, 2, &ok) ;      /* 2*nnz */
01081     c = COLAMD_C (n_col, &ok) ;     /* size of column structures */
01082     r = COLAMD_R (n_row, &ok) ;     /* size of row structures */
01083     s = t_add (s, c, &ok) ;
01084     s = t_add (s, r, &ok) ;
01085     s = t_add (s, n_col, &ok) ;     /* elbow room */
01086     s = t_add (s, nnz/5, &ok) ;     /* elbow room */
01087     ok = ok && (s < Int_MAX) ;
01088     return (ok ? s : 0) ;
01089 }
01090 
01091 
01092 /* ========================================================================== */
01093 /* === colamd_set_defaults ================================================== */
01094 /* ========================================================================== */
01095 
01096 /*
01097     The colamd_set_defaults routine sets the default values of the user-
01098     controllable parameters for colamd and symamd:
01099 
01100   Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
01101   entries are removed prior to ordering.  Columns with more than
01102   max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
01103   prior to ordering, and placed last in the output column ordering. 
01104 
01105   Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
01106   entries are removed prior to ordering, and placed last in the
01107   output ordering.
01108 
01109   knobs [0] dense row control
01110 
01111   knobs [1] dense column control
01112 
01113   knobs [2] if nonzero, do aggresive absorption
01114 
01115   knobs [3..19] unused, but future versions might use this
01116 
01117 */
01118 
01119 PUBLIC void COLAMD_set_defaults
01120 (
01121     /* === Parameters ======================================================= */
01122 
01123     double knobs [COLAMD_KNOBS]   /* knob array */
01124 )
01125 {
01126     /* === Local variables ================================================== */
01127 
01128     Int i ;
01129 
01130     if (!knobs)
01131     {
01132   return ;      /* no knobs to initialize */
01133     }
01134     for (i = 0 ; i < COLAMD_KNOBS ; i++)
01135     {
01136   knobs [i] = 0 ;
01137     }
01138     knobs [COLAMD_DENSE_ROW] = 10 ;
01139     knobs [COLAMD_DENSE_COL] = 10 ;
01140     knobs [COLAMD_AGGRESSIVE] = TRUE ;  /* default: do aggressive absorption*/
01141 }
01142 
01143 
01144 /* ========================================================================== */
01145 /* === symamd =============================================================== */
01146 /* ========================================================================== */
01147 
01148 PUBLIC Int SYMAMD_MAIN      /* return TRUE if OK, FALSE otherwise */
01149 (
01150     /* === Parameters ======================================================= */
01151 
01152     Int n,        /* number of rows and columns of A */
01153     Int A [],       /* row indices of A */
01154     Int p [],       /* column pointers of A */
01155     Int perm [],      /* output permutation, size n+1 */
01156     double knobs [COLAMD_KNOBS],  /* parameters (uses defaults if NULL) */
01157     Int stats [COLAMD_STATS],   /* output statistics and error codes */
01158     void * (*allocate) (size_t, size_t),
01159               /* pointer to calloc (ANSI C) or */
01160           /* mxCalloc (for MATLAB mexFunction) */
01161     void (*release) (void *)
01162               /* pointer to free (ANSI C) or */
01163               /* mxFree (for MATLAB mexFunction) */
01164 )
01165 {
01166     /* === Local variables ================================================== */
01167 
01168     Int *count ;    /* length of each column of M, and col pointer*/
01169     Int *mark ;     /* mark array for finding duplicate entries */
01170     Int *M ;      /* row indices of matrix M */
01171     size_t Mlen ;   /* length of M */
01172     Int n_row ;     /* number of rows in M */
01173     Int nnz ;     /* number of entries in A */
01174     Int i ;     /* row index of A */
01175     Int j ;     /* column index of A */
01176     Int k ;     /* row index of M */ 
01177     Int mnz ;     /* number of nonzeros in M */
01178     Int pp ;      /* index into a column of A */
01179     Int last_row ;    /* last row seen in the current column */
01180     Int length ;    /* number of nonzeros in a column */
01181 
01182     double cknobs [COLAMD_KNOBS] ;    /* knobs for colamd */
01183     double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */
01184 
01185 #ifndef NDEBUG
01186     colamd_get_debug ("symamd") ;
01187 #endif /* NDEBUG */
01188 
01189     /* === Check the input arguments ======================================== */
01190 
01191     if (!stats)
01192     {
01193   DEBUG0 (("symamd: stats not present\n")) ;
01194   return (FALSE) ;
01195     }
01196     for (i = 0 ; i < COLAMD_STATS ; i++)
01197     {
01198   stats [i] = 0 ;
01199     }
01200     stats [COLAMD_STATUS] = COLAMD_OK ;
01201     stats [COLAMD_INFO1] = -1 ;
01202     stats [COLAMD_INFO2] = -1 ;
01203 
01204     if (!A)
01205     {
01206       stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
01207   DEBUG0 (("symamd: A not present\n")) ;
01208   return (FALSE) ;
01209     }
01210 
01211     if (!p)   /* p is not present */
01212     {
01213   stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
01214   DEBUG0 (("symamd: p not present\n")) ;
01215       return (FALSE) ;
01216     }
01217 
01218     if (n < 0)    /* n must be >= 0 */
01219     {
01220   stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
01221   stats [COLAMD_INFO1] = n ;
01222   DEBUG0 (("symamd: n negative %d\n", n)) ;
01223       return (FALSE) ;
01224     }
01225 
01226     nnz = p [n] ;
01227     if (nnz < 0)  /* nnz must be >= 0 */
01228     {
01229   stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
01230   stats [COLAMD_INFO1] = nnz ;
01231   DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
01232   return (FALSE) ;
01233     }
01234 
01235     if (p [0] != 0)
01236     {
01237   stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
01238   stats [COLAMD_INFO1] = p [0] ;
01239   DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
01240   return (FALSE) ;
01241     }
01242 
01243     /* === If no knobs, set default knobs =================================== */
01244 
01245     if (!knobs)
01246     {
01247   COLAMD_set_defaults (default_knobs) ;
01248   knobs = default_knobs ;
01249     }
01250 
01251     /* === Allocate count and mark ========================================== */
01252 
01253     count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01254     if (!count)
01255     {
01256   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01257   DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
01258   return (FALSE) ;
01259     }
01260 
01261     mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01262     if (!mark)
01263     {
01264   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01265   (*release) ((void *) count) ;
01266   DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
01267   return (FALSE) ;
01268     }
01269 
01270     /* === Compute column counts of M, check if A is valid ================== */
01271 
01272     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
01273 
01274     for (i = 0 ; i < n ; i++)
01275     {
01276       mark [i] = -1 ;
01277     }
01278 
01279     for (j = 0 ; j < n ; j++)
01280     {
01281   last_row = -1 ;
01282 
01283   length = p [j+1] - p [j] ;
01284   if (length < 0)
01285   {
01286       /* column pointers must be non-decreasing */
01287       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
01288       stats [COLAMD_INFO1] = j ;
01289       stats [COLAMD_INFO2] = length ;
01290       (*release) ((void *) count) ;
01291       (*release) ((void *) mark) ;
01292       DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
01293       return (FALSE) ;
01294   }
01295 
01296   for (pp = p [j] ; pp < p [j+1] ; pp++)
01297   {
01298       i = A [pp] ;
01299       if (i < 0 || i >= n)
01300       {
01301     /* row index i, in column j, is out of bounds */
01302     stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
01303     stats [COLAMD_INFO1] = j ;
01304     stats [COLAMD_INFO2] = i ;
01305     stats [COLAMD_INFO3] = n ;
01306     (*release) ((void *) count) ;
01307     (*release) ((void *) mark) ;
01308     DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
01309     return (FALSE) ;
01310       }
01311 
01312       if (i <= last_row || mark [i] == j)
01313       {
01314     /* row index is unsorted or repeated (or both), thus col */
01315     /* is jumbled.  This is a notice, not an error condition. */
01316     stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
01317     stats [COLAMD_INFO1] = j ;
01318     stats [COLAMD_INFO2] = i ;
01319     (stats [COLAMD_INFO3]) ++ ;
01320     DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
01321       }
01322 
01323       if (i > j && mark [i] != j)
01324       {
01325     /* row k of M will contain column indices i and j */
01326     count [i]++ ;
01327     count [j]++ ;
01328       }
01329 
01330       /* mark the row as having been seen in this column */
01331       mark [i] = j ;
01332 
01333       last_row = i ;
01334   }
01335     }
01336 
01337     /* v2.4: removed free(mark) */
01338 
01339     /* === Compute column pointers of M ===================================== */
01340 
01341     /* use output permutation, perm, for column pointers of M */
01342     perm [0] = 0 ;
01343     for (j = 1 ; j <= n ; j++)
01344     {
01345   perm [j] = perm [j-1] + count [j-1] ;
01346     }
01347     for (j = 0 ; j < n ; j++)
01348     {
01349   count [j] = perm [j] ;
01350     }
01351 
01352     /* === Construct M ====================================================== */
01353 
01354     mnz = perm [n] ;
01355     n_row = mnz / 2 ;
01356     Mlen = COLAMD_recommended (mnz, n_row, n) ;
01357     M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
01358     DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
01359       n_row, n, mnz, (double) Mlen)) ;
01360 
01361     if (!M)
01362     {
01363   stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
01364   (*release) ((void *) count) ;
01365   (*release) ((void *) mark) ;
01366   DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
01367   return (FALSE) ;
01368     }
01369 
01370     k = 0 ;
01371 
01372     if (stats [COLAMD_STATUS] == COLAMD_OK)
01373     {
01374   /* Matrix is OK */
01375   for (j = 0 ; j < n ; j++)
01376   {
01377       ASSERT (p [j+1] - p [j] >= 0) ;
01378       for (pp = p [j] ; pp < p [j+1] ; pp++)
01379       {
01380     i = A [pp] ;
01381     ASSERT (i >= 0 && i < n) ;
01382     if (i > j)
01383     {
01384         /* row k of M contains column indices i and j */
01385         M [count [i]++] = k ;
01386         M [count [j]++] = k ;
01387         k++ ;
01388     }
01389       }
01390   }
01391     }
01392     else
01393     {
01394   /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
01395   DEBUG0 (("symamd: Duplicates in A.\n")) ;
01396   for (i = 0 ; i < n ; i++)
01397   {
01398       mark [i] = -1 ;
01399   }
01400   for (j = 0 ; j < n ; j++)
01401   {
01402       ASSERT (p [j+1] - p [j] >= 0) ;
01403       for (pp = p [j] ; pp < p [j+1] ; pp++)
01404       {
01405     i = A [pp] ;
01406     ASSERT (i >= 0 && i < n) ;
01407     if (i > j && mark [i] != j)
01408     {
01409         /* row k of M contains column indices i and j */
01410         M [count [i]++] = k ;
01411         M [count [j]++] = k ;
01412         k++ ;
01413         mark [i] = j ;
01414     }
01415       }
01416   }
01417   /* v2.4: free(mark) moved below */
01418     }
01419 
01420     /* count and mark no longer needed */
01421     (*release) ((void *) count) ;
01422     (*release) ((void *) mark) ;  /* v2.4: free (mark) moved here */
01423     ASSERT (k == n_row) ;
01424 
01425     /* === Adjust the knobs for M =========================================== */
01426 
01427     for (i = 0 ; i < COLAMD_KNOBS ; i++)
01428     {
01429   cknobs [i] = knobs [i] ;
01430     }
01431 
01432     /* there are no dense rows in M */
01433     cknobs [COLAMD_DENSE_ROW] = -1 ;
01434     cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
01435 
01436     /* === Order the columns of M =========================================== */
01437 
01438     /* v2.4: colamd cannot fail here, so the error check is removed */
01439     (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
01440 
01441     /* Note that the output permutation is now in perm */
01442 
01443     /* === get the statistics for symamd from colamd ======================== */
01444 
01445     /* a dense column in colamd means a dense row and col in symamd */
01446     stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
01447 
01448     /* === Free M =========================================================== */
01449 
01450     (*release) ((void *) M) ;
01451     DEBUG0 (("symamd: done.\n")) ;
01452     return (TRUE) ;
01453 
01454 }
01455 
01456 /* ========================================================================== */
01457 /* === colamd =============================================================== */
01458 /* ========================================================================== */
01459 
01460 /*
01461     The colamd routine computes a column ordering Q of a sparse matrix
01462     A such that the LU factorization P(AQ) = LU remains sparse, where P is
01463     selected via partial pivoting.   The routine can also be viewed as
01464     providing a permutation Q such that the Cholesky factorization
01465     (AQ)'(AQ) = LL' remains sparse.
01466 */
01467 
01468 PUBLIC Int COLAMD_MAIN    /* returns TRUE if successful, FALSE otherwise*/
01469 (
01470     /* === Parameters ======================================================= */
01471 
01472     Int n_row,      /* number of rows in A */
01473     Int n_col,      /* number of columns in A */
01474     Int Alen,     /* length of A */
01475     Int A [],     /* row indices of A */
01476     Int p [],     /* pointers to columns in A */
01477     double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
01478     Int stats [COLAMD_STATS]  /* output statistics and error codes */
01479 )
01480 {
01481     /* === Local variables ================================================== */
01482 
01483     Int i ;     /* loop index */
01484     Int nnz ;     /* nonzeros in A */
01485     size_t Row_size ;   /* size of Row [], in integers */
01486     size_t Col_size ;   /* size of Col [], in integers */
01487     size_t need ;   /* minimum required length of A */
01488     Colamd_Row *Row ;   /* pointer into A of Row [0..n_row] array */
01489     Colamd_Col *Col ;   /* pointer into A of Col [0..n_col] array */
01490     Int n_col2 ;    /* number of non-dense, non-empty columns */
01491     Int n_row2 ;    /* number of non-dense, non-empty rows */
01492     Int ngarbage ;    /* number of garbage collections performed */
01493     Int max_deg ;   /* maximum row degree */
01494     double default_knobs [COLAMD_KNOBS] ; /* default knobs array */
01495     Int aggressive ;    /* do aggressive absorption */
01496     int ok ;
01497 
01498 #ifndef NDEBUG
01499     colamd_get_debug ("colamd") ;
01500 #endif /* NDEBUG */
01501 
01502     /* === Check the input arguments ======================================== */
01503 
01504     if (!stats)
01505     {
01506   DEBUG0 (("colamd: stats not present\n")) ;
01507   return (FALSE) ;
01508     }
01509     for (i = 0 ; i < COLAMD_STATS ; i++)
01510     {
01511   stats [i] = 0 ;
01512     }
01513     stats [COLAMD_STATUS] = COLAMD_OK ;
01514     stats [COLAMD_INFO1] = -1 ;
01515     stats [COLAMD_INFO2] = -1 ;
01516 
01517     if (!A)   /* A is not present */
01518     {
01519   stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
01520   DEBUG0 (("colamd: A not present\n")) ;
01521   return (FALSE) ;
01522     }
01523 
01524     if (!p)   /* p is not present */
01525     {
01526   stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
01527   DEBUG0 (("colamd: p not present\n")) ;
01528       return (FALSE) ;
01529     }
01530 
01531     if (n_row < 0)  /* n_row must be >= 0 */
01532     {
01533   stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
01534   stats [COLAMD_INFO1] = n_row ;
01535   DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
01536       return (FALSE) ;
01537     }
01538 
01539     if (n_col < 0)  /* n_col must be >= 0 */
01540     {
01541   stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
01542   stats [COLAMD_INFO1] = n_col ;
01543   DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
01544       return (FALSE) ;
01545     }
01546 
01547     nnz = p [n_col] ;
01548     if (nnz < 0)  /* nnz must be >= 0 */
01549     {
01550   stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
01551   stats [COLAMD_INFO1] = nnz ;
01552   DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
01553   return (FALSE) ;
01554     }
01555 
01556     if (p [0] != 0)
01557     {
01558   stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
01559   stats [COLAMD_INFO1] = p [0] ;
01560   DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
01561   return (FALSE) ;
01562     }
01563 
01564     /* === If no knobs, set default knobs =================================== */
01565 
01566     if (!knobs)
01567     {
01568   COLAMD_set_defaults (default_knobs) ;
01569   knobs = default_knobs ;
01570     }
01571 
01572     aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
01573 
01574     /* === Allocate the Row and Col arrays from array A ===================== */
01575 
01576     ok = TRUE ;
01577     Col_size = COLAMD_C (n_col, &ok) ;      /* size of Col array of structs */
01578     Row_size = COLAMD_R (n_row, &ok) ;      /* size of Row array of structs */
01579 
01580     /* need = 2*nnz + n_col + Col_size + Row_size ; */
01581     need = t_mult (nnz, 2, &ok) ;
01582     need = t_add (need, n_col, &ok) ;
01583     need = t_add (need, Col_size, &ok) ;
01584     need = t_add (need, Row_size, &ok) ;
01585 
01586     if (!ok || need > (size_t) Alen || need > Int_MAX)
01587     {
01588   /* not enough space in array A to perform the ordering */
01589   stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
01590   stats [COLAMD_INFO1] = need ;
01591   stats [COLAMD_INFO2] = Alen ;
01592   DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
01593   return (FALSE) ;
01594     }
01595 
01596     Alen -= Col_size + Row_size ;
01597     Col = (Colamd_Col *) &A [Alen] ;
01598     Row = (Colamd_Row *) &A [Alen + Col_size] ;
01599 
01600     /* === Construct the row and column data structures ===================== */
01601 
01602     if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
01603     {
01604   /* input matrix is invalid */
01605   DEBUG0 (("colamd: Matrix invalid\n")) ;
01606   return (FALSE) ;
01607     }
01608 
01609     /* === Initialize scores, kill dense rows/columns ======================= */
01610 
01611     init_scoring (n_row, n_col, Row, Col, A, p, knobs,
01612   &n_row2, &n_col2, &max_deg) ;
01613 
01614     /* === Order the supercolumns =========================================== */
01615 
01616     ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
01617   n_col2, max_deg, 2*nnz, aggressive) ;
01618 
01619     /* === Order the non-principal columns ================================== */
01620 
01621     order_children (n_col, Col, p) ;
01622 
01623     /* === Return statistics in stats ======================================= */
01624 
01625     stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
01626     stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
01627     stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
01628     DEBUG0 (("colamd: done.\n")) ; 
01629     return (TRUE) ;
01630 }
01631 
01632 
01633 /* ========================================================================== */
01634 /* === colamd_report ======================================================== */
01635 /* ========================================================================== */
01636 
01637 PUBLIC void COLAMD_report
01638 (
01639     Int stats [COLAMD_STATS]
01640 )
01641 {
01642     print_report ("colamd", stats) ;
01643 }
01644 
01645 
01646 /* ========================================================================== */
01647 /* === symamd_report ======================================================== */
01648 /* ========================================================================== */
01649 
01650 PUBLIC void SYMAMD_report
01651 (
01652     Int stats [COLAMD_STATS]
01653 )
01654 {
01655     print_report ("symamd", stats) ;
01656 }
01657 
01658 
01659 
01660 /* ========================================================================== */
01661 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
01662 /* ========================================================================== */
01663 
01664 /* There are no user-callable routines beyond this point in the file */
01665 
01666 
01667 /* ========================================================================== */
01668 /* === init_rows_cols ======================================================= */
01669 /* ========================================================================== */
01670 
01671 /*
01672     Takes the column form of the matrix in A and creates the row form of the
01673     matrix.  Also, row and column attributes are stored in the Col and Row
01674     structs.  If the columns are un-sorted or contain duplicate row indices,
01675     this routine will also sort and remove duplicate row indices from the
01676     column form of the matrix.  Returns FALSE if the matrix is invalid,
01677     TRUE otherwise.  Not user-callable.
01678 */
01679 
01680 PRIVATE Int init_rows_cols  /* returns TRUE if OK, or FALSE otherwise */
01681 (
01682     /* === Parameters ======================================================= */
01683 
01684     Int n_row,      /* number of rows of A */
01685     Int n_col,      /* number of columns of A */
01686     Colamd_Row Row [],    /* of size n_row+1 */
01687     Colamd_Col Col [],    /* of size n_col+1 */
01688     Int A [],     /* row indices of A, of size Alen */
01689     Int p [],     /* pointers to columns in A, of size n_col+1 */
01690     Int stats [COLAMD_STATS]  /* colamd statistics */ 
01691 )
01692 {
01693     /* === Local variables ================================================== */
01694 
01695     Int col ;     /* a column index */
01696     Int row ;     /* a row index */
01697     Int *cp ;     /* a column pointer */
01698     Int *cp_end ;   /* a pointer to the end of a column */
01699     Int *rp ;     /* a row pointer */
01700     Int *rp_end ;   /* a pointer to the end of a row */
01701     Int last_row ;    /* previous row */
01702 
01703     /* === Initialize columns, and check column pointers ==================== */
01704 
01705     for (col = 0 ; col < n_col ; col++)
01706     {
01707   Col [col].start = p [col] ;
01708   Col [col].length = p [col+1] - p [col] ;
01709 
01710   if (Col [col].length < 0)
01711   {
01712       /* column pointers must be non-decreasing */
01713       stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
01714       stats [COLAMD_INFO1] = col ;
01715       stats [COLAMD_INFO2] = Col [col].length ;
01716       DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
01717       return (FALSE) ;
01718   }
01719 
01720   Col [col].shared1.thickness = 1 ;
01721   Col [col].shared2.score = 0 ;
01722   Col [col].shared3.prev = EMPTY ;
01723   Col [col].shared4.degree_next = EMPTY ;
01724     }
01725 
01726     /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
01727 
01728     /* === Scan columns, compute row degrees, and check row indices ========= */
01729 
01730     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
01731 
01732     for (row = 0 ; row < n_row ; row++)
01733     {
01734   Row [row].length = 0 ;
01735   Row [row].shared2.mark = -1 ;
01736     }
01737 
01738     for (col = 0 ; col < n_col ; col++)
01739     {
01740   last_row = -1 ;
01741 
01742   cp = &A [p [col]] ;
01743   cp_end = &A [p [col+1]] ;
01744 
01745   while (cp < cp_end)
01746   {
01747       row = *cp++ ;
01748 
01749       /* make sure row indices within range */
01750       if (row < 0 || row >= n_row)
01751       {
01752     stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
01753     stats [COLAMD_INFO1] = col ;
01754     stats [COLAMD_INFO2] = row ;
01755     stats [COLAMD_INFO3] = n_row ;
01756     DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
01757     return (FALSE) ;
01758       }
01759 
01760       if (row <= last_row || Row [row].shared2.mark == col)
01761       {
01762     /* row index are unsorted or repeated (or both), thus col */
01763     /* is jumbled.  This is a notice, not an error condition. */
01764     stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
01765     stats [COLAMD_INFO1] = col ;
01766     stats [COLAMD_INFO2] = row ;
01767     (stats [COLAMD_INFO3]) ++ ;
01768     DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
01769       }
01770 
01771       if (Row [row].shared2.mark != col)
01772       {
01773     Row [row].length++ ;
01774       }
01775       else
01776       {
01777     /* this is a repeated entry in the column, */
01778     /* it will be removed */
01779     Col [col].length-- ;
01780       }
01781 
01782       /* mark the row as having been seen in this column */
01783       Row [row].shared2.mark = col ;
01784 
01785       last_row = row ;
01786   }
01787     }
01788 
01789     /* === Compute row pointers ============================================= */
01790 
01791     /* row form of the matrix starts directly after the column */
01792     /* form of matrix in A */
01793     Row [0].start = p [n_col] ;
01794     Row [0].shared1.p = Row [0].start ;
01795     Row [0].shared2.mark = -1 ;
01796     for (row = 1 ; row < n_row ; row++)
01797     {
01798   Row [row].start = Row [row-1].start + Row [row-1].length ;
01799   Row [row].shared1.p = Row [row].start ;
01800   Row [row].shared2.mark = -1 ;
01801     }
01802 
01803     /* === Create row form ================================================== */
01804 
01805     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
01806     {
01807   /* if cols jumbled, watch for repeated row indices */
01808   for (col = 0 ; col < n_col ; col++)
01809   {
01810       cp = &A [p [col]] ;
01811       cp_end = &A [p [col+1]] ;
01812       while (cp < cp_end)
01813       {
01814     row = *cp++ ;
01815     if (Row [row].shared2.mark != col)
01816     {
01817         A [(Row [row].shared1.p)++] = col ;
01818         Row [row].shared2.mark = col ;
01819     }
01820       }
01821   }
01822     }
01823     else
01824     {
01825   /* if cols not jumbled, we don't need the mark (this is faster) */
01826   for (col = 0 ; col < n_col ; col++)
01827   {
01828       cp = &A [p [col]] ;
01829       cp_end = &A [p [col+1]] ;
01830       while (cp < cp_end)
01831       {
01832     A [(Row [*cp++].shared1.p)++] = col ;
01833       }
01834   }
01835     }
01836 
01837     /* === Clear the row marks and set row degrees ========================== */
01838 
01839     for (row = 0 ; row < n_row ; row++)
01840     {
01841   Row [row].shared2.mark = 0 ;
01842   Row [row].shared1.degree = Row [row].length ;
01843     }
01844 
01845     /* === See if we need to re-create columns ============================== */
01846 
01847     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
01848     {
01849       DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
01850 
01851 #ifndef NDEBUG
01852   /* make sure column lengths are correct */
01853   for (col = 0 ; col < n_col ; col++)
01854   {
01855       p [col] = Col [col].length ;
01856   }
01857   for (row = 0 ; row < n_row ; row++)
01858   {
01859       rp = &A [Row [row].start] ;
01860       rp_end = rp + Row [row].length ;
01861       while (rp < rp_end)
01862       {
01863     p [*rp++]-- ;
01864       }
01865   }
01866   for (col = 0 ; col < n_col ; col++)
01867   {
01868       ASSERT (p [col] == 0) ;
01869   }
01870   /* now p is all zero (different than when debugging is turned off) */
01871 #endif /* NDEBUG */
01872 
01873   /* === Compute col pointers ========================================= */
01874 
01875   /* col form of the matrix starts at A [0]. */
01876   /* Note, we may have a gap between the col form and the row */
01877   /* form if there were duplicate entries, if so, it will be */
01878   /* removed upon the first garbage collection */
01879   Col [0].start = 0 ;
01880   p [0] = Col [0].start ;
01881   for (col = 1 ; col < n_col ; col++)
01882   {
01883       /* note that the lengths here are for pruned columns, i.e. */
01884       /* no duplicate row indices will exist for these columns */
01885       Col [col].start = Col [col-1].start + Col [col-1].length ;
01886       p [col] = Col [col].start ;
01887   }
01888 
01889   /* === Re-create col form =========================================== */
01890 
01891   for (row = 0 ; row < n_row ; row++)
01892   {
01893       rp = &A [Row [row].start] ;
01894       rp_end = rp + Row [row].length ;
01895       while (rp < rp_end)
01896       {
01897     A [(p [*rp++])++] = row ;
01898       }
01899   }
01900     }
01901 
01902     /* === Done.  Matrix is not (or no longer) jumbled ====================== */
01903 
01904     return (TRUE) ;
01905 }
01906 
01907 
01908 /* ========================================================================== */
01909 /* === init_scoring ========================================================= */
01910 /* ========================================================================== */
01911 
01912 /*
01913     Kills dense or empty columns and rows, calculates an initial score for
01914     each column, and places all columns in the degree lists.  Not user-callable.
01915 */
01916 
01917 PRIVATE void init_scoring
01918 (
01919     /* === Parameters ======================================================= */
01920 
01921     Int n_row,      /* number of rows of A */
01922     Int n_col,      /* number of columns of A */
01923     Colamd_Row Row [],    /* of size n_row+1 */
01924     Colamd_Col Col [],    /* of size n_col+1 */
01925     Int A [],     /* column form and row form of A */
01926     Int head [],    /* of size n_col+1 */
01927     double knobs [COLAMD_KNOBS],/* parameters */
01928     Int *p_n_row2,    /* number of non-dense, non-empty rows */
01929     Int *p_n_col2,    /* number of non-dense, non-empty columns */
01930     Int *p_max_deg    /* maximum row degree */
01931 )
01932 {
01933     /* === Local variables ================================================== */
01934 
01935     Int c ;     /* a column index */
01936     Int r, row ;    /* a row index */
01937     Int *cp ;     /* a column pointer */
01938     Int deg ;     /* degree of a row or column */
01939     Int *cp_end ;   /* a pointer to the end of a column */
01940     Int *new_cp ;   /* new column pointer */
01941     Int col_length ;    /* length of pruned column */
01942     Int score ;     /* current column score */
01943     Int n_col2 ;    /* number of non-dense, non-empty columns */
01944     Int n_row2 ;    /* number of non-dense, non-empty rows */
01945     Int dense_row_count ; /* remove rows with more entries than this */
01946     Int dense_col_count ; /* remove cols with more entries than this */
01947     Int min_score ;   /* smallest column score */
01948     Int max_deg ;   /* maximum row degree */
01949     Int next_col ;    /* Used to add to degree list.*/
01950 
01951 #ifndef NDEBUG
01952     Int debug_count ;   /* debug only. */
01953 #endif /* NDEBUG */
01954 
01955     /* === Extract knobs ==================================================== */
01956 
01957     /* Note: if knobs contains a NaN, this is undefined: */
01958     if (knobs [COLAMD_DENSE_ROW] < 0)
01959     {
01960   /* only remove completely dense rows */
01961   dense_row_count = n_col-1 ;
01962     }
01963     else
01964     {
01965   dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
01966     }
01967     if (knobs [COLAMD_DENSE_COL] < 0)
01968     {
01969   /* only remove completely dense columns */
01970   dense_col_count = n_row-1 ;
01971     }
01972     else
01973     {
01974   dense_col_count =
01975       DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
01976     }
01977 
01978     DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
01979     max_deg = 0 ;
01980     n_col2 = n_col ;
01981     n_row2 = n_row ;
01982 
01983     /* === Kill empty columns =============================================== */
01984 
01985     /* Put the empty columns at the end in their natural order, so that LU */
01986     /* factorization can proceed as far as possible. */
01987     for (c = n_col-1 ; c >= 0 ; c--)
01988     {
01989   deg = Col [c].length ;
01990   if (deg == 0)
01991   {
01992       /* this is a empty column, kill and order it last */
01993       Col [c].shared2.order = --n_col2 ;
01994       KILL_PRINCIPAL_COL (c) ;
01995   }
01996     }
01997     DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
01998 
01999     /* === Kill dense columns =============================================== */
02000 
02001     /* Put the dense columns at the end, in their natural order */
02002     for (c = n_col-1 ; c >= 0 ; c--)
02003     {
02004   /* skip any dead columns */
02005   if (COL_IS_DEAD (c))
02006   {
02007       continue ;
02008   }
02009   deg = Col [c].length ;
02010   if (deg > dense_col_count)
02011   {
02012       /* this is a dense column, kill and order it last */
02013       Col [c].shared2.order = --n_col2 ;
02014       /* decrement the row degrees */
02015       cp = &A [Col [c].start] ;
02016       cp_end = cp + Col [c].length ;
02017       while (cp < cp_end)
02018       {
02019     Row [*cp++].shared1.degree-- ;
02020       }
02021       KILL_PRINCIPAL_COL (c) ;
02022   }
02023     }
02024     DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
02025 
02026     /* === Kill dense and empty rows ======================================== */
02027 
02028     for (r = 0 ; r < n_row ; r++)
02029     {
02030   deg = Row [r].shared1.degree ;
02031   ASSERT (deg >= 0 && deg <= n_col) ;
02032   if (deg > dense_row_count || deg == 0)
02033   {
02034       /* kill a dense or empty row */
02035       KILL_ROW (r) ;
02036       --n_row2 ;
02037   }
02038   else
02039   {
02040       /* keep track of max degree of remaining rows */
02041       max_deg = MAX (max_deg, deg) ;
02042   }
02043     }
02044     DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
02045 
02046     /* === Compute initial column scores ==================================== */
02047 
02048     /* At this point the row degrees are accurate.  They reflect the number */
02049     /* of "live" (non-dense) columns in each row.  No empty rows exist. */
02050     /* Some "live" columns may contain only dead rows, however.  These are */
02051     /* pruned in the code below. */
02052 
02053     /* now find the initial matlab score for each column */
02054     for (c = n_col-1 ; c >= 0 ; c--)
02055     {
02056   /* skip dead column */
02057   if (COL_IS_DEAD (c))
02058   {
02059       continue ;
02060   }
02061   score = 0 ;
02062   cp = &A [Col [c].start] ;
02063   new_cp = cp ;
02064   cp_end = cp + Col [c].length ;
02065   while (cp < cp_end)
02066   {
02067       /* get a row */
02068       row = *cp++ ;
02069       /* skip if dead */
02070       if (ROW_IS_DEAD (row))
02071       {
02072     continue ;
02073       }
02074       /* compact the column */
02075       *new_cp++ = row ;
02076       /* add row's external degree */
02077       score += Row [row].shared1.degree - 1 ;
02078       /* guard against integer overflow */
02079       score = MIN (score, n_col) ;
02080   }
02081   /* determine pruned column length */
02082   col_length = (Int) (new_cp - &A [Col [c].start]) ;
02083   if (col_length == 0)
02084   {
02085       /* a newly-made null column (all rows in this col are "dense" */
02086       /* and have already been killed) */
02087       DEBUG2 (("Newly null killed: %d\n", c)) ;
02088       Col [c].shared2.order = --n_col2 ;
02089       KILL_PRINCIPAL_COL (c) ;
02090   }
02091   else
02092   {
02093       /* set column length and set score */
02094       ASSERT (score >= 0) ;
02095       ASSERT (score <= n_col) ;
02096       Col [c].length = col_length ;
02097       Col [c].shared2.score = score ;
02098   }
02099     }
02100     DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
02101       n_col-n_col2)) ;
02102 
02103     /* At this point, all empty rows and columns are dead.  All live columns */
02104     /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
02105     /* yet).  Rows may contain dead columns, but all live rows contain at */
02106     /* least one live column. */
02107 
02108 #ifndef NDEBUG
02109     debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
02110 #endif /* NDEBUG */
02111 
02112     /* === Initialize degree lists ========================================== */
02113 
02114 #ifndef NDEBUG
02115     debug_count = 0 ;
02116 #endif /* NDEBUG */
02117 
02118     /* clear the hash buckets */
02119     for (c = 0 ; c <= n_col ; c++)
02120     {
02121   head [c] = EMPTY ;
02122     }
02123     min_score = n_col ;
02124     /* place in reverse order, so low column indices are at the front */
02125     /* of the lists.  This is to encourage natural tie-breaking */
02126     for (c = n_col-1 ; c >= 0 ; c--)
02127     {
02128   /* only add principal columns to degree lists */
02129   if (COL_IS_ALIVE (c))
02130   {
02131       DEBUG4 (("place %d score %d minscore %d ncol %d\n",
02132     c, Col [c].shared2.score, min_score, n_col)) ;
02133 
02134       /* === Add columns score to DList =============================== */
02135 
02136       score = Col [c].shared2.score ;
02137 
02138       ASSERT (min_score >= 0) ;
02139       ASSERT (min_score <= n_col) ;
02140       ASSERT (score >= 0) ;
02141       ASSERT (score <= n_col) ;
02142       ASSERT (head [score] >= EMPTY) ;
02143 
02144       /* now add this column to dList at proper score location */
02145       next_col = head [score] ;
02146       Col [c].shared3.prev = EMPTY ;
02147       Col [c].shared4.degree_next = next_col ;
02148 
02149       /* if there already was a column with the same score, set its */
02150       /* previous pointer to this new column */
02151       if (next_col != EMPTY)
02152       {
02153     Col [next_col].shared3.prev = c ;
02154       }
02155       head [score] = c ;
02156 
02157       /* see if this score is less than current min */
02158       min_score = MIN (min_score, score) ;
02159 
02160 #ifndef NDEBUG
02161       debug_count++ ;
02162 #endif /* NDEBUG */
02163 
02164   }
02165     }
02166 
02167 #ifndef NDEBUG
02168     DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
02169   debug_count, n_col, n_col-debug_count)) ;
02170     ASSERT (debug_count == n_col2) ;
02171     debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
02172 #endif /* NDEBUG */
02173 
02174     /* === Return number of remaining columns, and max row degree =========== */
02175 
02176     *p_n_col2 = n_col2 ;
02177     *p_n_row2 = n_row2 ;
02178     *p_max_deg = max_deg ;
02179 }
02180 
02181 
02182 /* ========================================================================== */
02183 /* === find_ordering ======================================================== */
02184 /* ========================================================================== */
02185 
02186 /*
02187     Order the principal columns of the supercolumn form of the matrix
02188     (no supercolumns on input).  Uses a minimum approximate column minimum
02189     degree ordering method.  Not user-callable.
02190 */
02191 
02192 PRIVATE Int find_ordering /* return the number of garbage collections */
02193 (
02194     /* === Parameters ======================================================= */
02195 
02196     Int n_row,      /* number of rows of A */
02197     Int n_col,      /* number of columns of A */
02198     Int Alen,     /* size of A, 2*nnz + n_col or larger */
02199     Colamd_Row Row [],    /* of size n_row+1 */
02200     Colamd_Col Col [],    /* of size n_col+1 */
02201     Int A [],     /* column form and row form of A */
02202     Int head [],    /* of size n_col+1 */
02203     Int n_col2,     /* Remaining columns to order */
02204     Int max_deg,    /* Maximum row degree */
02205     Int pfree,      /* index of first free slot (2*nnz on entry) */
02206     Int aggressive
02207 )
02208 {
02209     /* === Local variables ================================================== */
02210 
02211     Int k ;     /* current pivot ordering step */
02212     Int pivot_col ;   /* current pivot column */
02213     Int *cp ;     /* a column pointer */
02214     Int *rp ;     /* a row pointer */
02215     Int pivot_row ;   /* current pivot row */
02216     Int *new_cp ;   /* modified column pointer */
02217     Int *new_rp ;   /* modified row pointer */
02218     Int pivot_row_start ; /* pointer to start of pivot row */
02219     Int pivot_row_degree ;  /* number of columns in pivot row */
02220     Int pivot_row_length ;  /* number of supercolumns in pivot row */
02221     Int pivot_col_score ; /* score of pivot column */
02222     Int needed_memory ;   /* free space needed for pivot row */
02223     Int *cp_end ;   /* pointer to the end of a column */
02224     Int *rp_end ;   /* pointer to the end of a row */
02225     Int row ;     /* a row index */
02226     Int col ;     /* a column index */
02227     Int max_score ;   /* maximum possible score */
02228     Int cur_score ;   /* score of current column */
02229     unsigned Int hash ;   /* hash value for supernode detection */
02230     Int head_column ;   /* head of hash bucket */
02231     Int first_col ;   /* first column in hash bucket */
02232     Int tag_mark ;    /* marker value for mark array */
02233     Int row_mark ;    /* Row [row].shared2.mark */
02234     Int set_difference ;  /* set difference size of row with pivot row */
02235     Int min_score ;   /* smallest column score */
02236     Int col_thickness ;   /* "thickness" (no. of columns in a supercol) */
02237     Int max_mark ;    /* maximum value of tag_mark */
02238     Int pivot_col_thickness ; /* number of columns represented by pivot col */
02239     Int prev_col ;    /* Used by Dlist operations. */
02240     Int next_col ;    /* Used by Dlist operations. */
02241     Int ngarbage ;    /* number of garbage collections performed */
02242 
02243 #ifndef NDEBUG
02244     Int debug_d ;   /* debug loop counter */
02245     Int debug_step = 0 ;  /* debug loop counter */
02246 #endif /* NDEBUG */
02247 
02248     /* === Initialization and clear mark ==================================== */
02249 
02250     max_mark = INT_MAX - n_col ;  /* INT_MAX defined in <limits.h> */
02251     tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02252     min_score = 0 ;
02253     ngarbage = 0 ;
02254     DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
02255 
02256     /* === Order the columns ================================================ */
02257 
02258     for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
02259     {
02260 
02261 #ifndef NDEBUG
02262   if (debug_step % 100 == 0)
02263   {
02264       DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
02265   }
02266   else
02267   {
02268       DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
02269   }
02270   debug_step++ ;
02271   debug_deg_lists (n_row, n_col, Row, Col, head,
02272     min_score, n_col2-k, max_deg) ;
02273   debug_matrix (n_row, n_col, Row, Col, A) ;
02274 #endif /* NDEBUG */
02275 
02276   /* === Select pivot column, and order it ============================ */
02277 
02278   /* make sure degree list isn't empty */
02279   ASSERT (min_score >= 0) ;
02280   ASSERT (min_score <= n_col) ;
02281   ASSERT (head [min_score] >= EMPTY) ;
02282 
02283 #ifndef NDEBUG
02284   for (debug_d = 0 ; debug_d < min_score ; debug_d++)
02285   {
02286       ASSERT (head [debug_d] == EMPTY) ;
02287   }
02288 #endif /* NDEBUG */
02289 
02290   /* get pivot column from head of minimum degree list */
02291   while (head [min_score] == EMPTY && min_score < n_col)
02292   {
02293       min_score++ ;
02294   }
02295   pivot_col = head [min_score] ;
02296   ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
02297   next_col = Col [pivot_col].shared4.degree_next ;
02298   head [min_score] = next_col ;
02299   if (next_col != EMPTY)
02300   {
02301       Col [next_col].shared3.prev = EMPTY ;
02302   }
02303 
02304   ASSERT (COL_IS_ALIVE (pivot_col)) ;
02305 
02306   /* remember score for defrag check */
02307   pivot_col_score = Col [pivot_col].shared2.score ;
02308 
02309   /* the pivot column is the kth column in the pivot order */
02310   Col [pivot_col].shared2.order = k ;
02311 
02312   /* increment order count by column thickness */
02313   pivot_col_thickness = Col [pivot_col].shared1.thickness ;
02314   k += pivot_col_thickness ;
02315   ASSERT (pivot_col_thickness > 0) ;
02316   DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
02317 
02318   /* === Garbage_collection, if necessary ============================= */
02319 
02320   needed_memory = MIN (pivot_col_score, n_col - k) ;
02321   if (pfree + needed_memory >= Alen)
02322   {
02323       pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
02324       ngarbage++ ;
02325       /* after garbage collection we will have enough */
02326       ASSERT (pfree + needed_memory < Alen) ;
02327       /* garbage collection has wiped out the Row[].shared2.mark array */
02328       tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02329 
02330 #ifndef NDEBUG
02331       debug_matrix (n_row, n_col, Row, Col, A) ;
02332 #endif /* NDEBUG */
02333   }
02334 
02335   /* === Compute pivot row pattern ==================================== */
02336 
02337   /* get starting location for this new merged row */
02338   pivot_row_start = pfree ;
02339 
02340   /* initialize new row counts to zero */
02341   pivot_row_degree = 0 ;
02342 
02343   /* tag pivot column as having been visited so it isn't included */
02344   /* in merged pivot row */
02345   Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
02346 
02347   /* pivot row is the union of all rows in the pivot column pattern */
02348   cp = &A [Col [pivot_col].start] ;
02349   cp_end = cp + Col [pivot_col].length ;
02350   while (cp < cp_end)
02351   {
02352       /* get a row */
02353       row = *cp++ ;
02354       DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
02355       /* skip if row is dead */
02356       if (ROW_IS_ALIVE (row))
02357       {
02358     rp = &A [Row [row].start] ;
02359     rp_end = rp + Row [row].length ;
02360     while (rp < rp_end)
02361     {
02362         /* get a column */
02363         col = *rp++ ;
02364         /* add the column, if alive and untagged */
02365         col_thickness = Col [col].shared1.thickness ;
02366         if (col_thickness > 0 && COL_IS_ALIVE (col))
02367         {
02368       /* tag column in pivot row */
02369       Col [col].shared1.thickness = -col_thickness ;
02370       ASSERT (pfree < Alen) ;
02371       /* place column in pivot row */
02372       A [pfree++] = col ;
02373       pivot_row_degree += col_thickness ;
02374         }
02375     }
02376       }
02377   }
02378 
02379   /* clear tag on pivot column */
02380   Col [pivot_col].shared1.thickness = pivot_col_thickness ;
02381   max_deg = MAX (max_deg, pivot_row_degree) ;
02382 
02383 #ifndef NDEBUG
02384   DEBUG3 (("check2\n")) ;
02385   debug_mark (n_row, Row, tag_mark, max_mark) ;
02386 #endif /* NDEBUG */
02387 
02388   /* === Kill all rows used to construct pivot row ==================== */
02389 
02390   /* also kill pivot row, temporarily */
02391   cp = &A [Col [pivot_col].start] ;
02392   cp_end = cp + Col [pivot_col].length ;
02393   while (cp < cp_end)
02394   {
02395       /* may be killing an already dead row */
02396       row = *cp++ ;
02397       DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
02398       KILL_ROW (row) ;
02399   }
02400 
02401   /* === Select a row index to use as the new pivot row =============== */
02402 
02403   pivot_row_length = pfree - pivot_row_start ;
02404   if (pivot_row_length > 0)
02405   {
02406       /* pick the "pivot" row arbitrarily (first row in col) */
02407       pivot_row = A [Col [pivot_col].start] ;
02408       DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
02409   }
02410   else
02411   {
02412       /* there is no pivot row, since it is of zero length */
02413       pivot_row = EMPTY ;
02414       ASSERT (pivot_row_length == 0) ;
02415   }
02416   ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
02417 
02418   /* === Approximate degree computation =============================== */
02419 
02420   /* Here begins the computation of the approximate degree.  The column */
02421   /* score is the sum of the pivot row "length", plus the size of the */
02422   /* set differences of each row in the column minus the pattern of the */
02423   /* pivot row itself.  The column ("thickness") itself is also */
02424   /* excluded from the column score (we thus use an approximate */
02425   /* external degree). */
02426 
02427   /* The time taken by the following code (compute set differences, and */
02428   /* add them up) is proportional to the size of the data structure */
02429   /* being scanned - that is, the sum of the sizes of each column in */
02430   /* the pivot row.  Thus, the amortized time to compute a column score */
02431   /* is proportional to the size of that column (where size, in this */
02432   /* context, is the column "length", or the number of row indices */
02433   /* in that column).  The number of row indices in a column is */
02434   /* monotonically non-decreasing, from the length of the original */
02435   /* column on input to colamd. */
02436 
02437   /* === Compute set differences ====================================== */
02438 
02439   DEBUG3 (("** Computing set differences phase. **\n")) ;
02440 
02441   /* pivot row is currently dead - it will be revived later. */
02442 
02443   DEBUG3 (("Pivot row: ")) ;
02444   /* for each column in pivot row */
02445   rp = &A [pivot_row_start] ;
02446   rp_end = rp + pivot_row_length ;
02447   while (rp < rp_end)
02448   {
02449       col = *rp++ ;
02450       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
02451       DEBUG3 (("Col: %d\n", col)) ;
02452 
02453       /* clear tags used to construct pivot row pattern */
02454       col_thickness = -Col [col].shared1.thickness ;
02455       ASSERT (col_thickness > 0) ;
02456       Col [col].shared1.thickness = col_thickness ;
02457 
02458       /* === Remove column from degree list =========================== */
02459 
02460       cur_score = Col [col].shared2.score ;
02461       prev_col = Col [col].shared3.prev ;
02462       next_col = Col [col].shared4.degree_next ;
02463       ASSERT (cur_score >= 0) ;
02464       ASSERT (cur_score <= n_col) ;
02465       ASSERT (cur_score >= EMPTY) ;
02466       if (prev_col == EMPTY)
02467       {
02468     head [cur_score] = next_col ;
02469       }
02470       else
02471       {
02472     Col [prev_col].shared4.degree_next = next_col ;
02473       }
02474       if (next_col != EMPTY)
02475       {
02476     Col [next_col].shared3.prev = prev_col ;
02477       }
02478 
02479       /* === Scan the column ========================================== */
02480 
02481       cp = &A [Col [col].start] ;
02482       cp_end = cp + Col [col].length ;
02483       while (cp < cp_end)
02484       {
02485     /* get a row */
02486     row = *cp++ ;
02487     row_mark = Row [row].shared2.mark ;
02488     /* skip if dead */
02489     if (ROW_IS_MARKED_DEAD (row_mark))
02490     {
02491         continue ;
02492     }
02493     ASSERT (row != pivot_row) ;
02494     set_difference = row_mark - tag_mark ;
02495     /* check if the row has been seen yet */
02496     if (set_difference < 0)
02497     {
02498         ASSERT (Row [row].shared1.degree <= max_deg) ;
02499         set_difference = Row [row].shared1.degree ;
02500     }
02501     /* subtract column thickness from this row's set difference */
02502     set_difference -= col_thickness ;
02503     ASSERT (set_difference >= 0) ;
02504     /* absorb this row if the set difference becomes zero */
02505     if (set_difference == 0 && aggressive)
02506     {
02507         DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
02508         KILL_ROW (row) ;
02509     }
02510     else
02511     {
02512         /* save the new mark */
02513         Row [row].shared2.mark = set_difference + tag_mark ;
02514     }
02515       }
02516   }
02517 
02518 #ifndef NDEBUG
02519   debug_deg_lists (n_row, n_col, Row, Col, head,
02520     min_score, n_col2-k-pivot_row_degree, max_deg) ;
02521 #endif /* NDEBUG */
02522 
02523   /* === Add up set differences for each column ======================= */
02524 
02525   DEBUG3 (("** Adding set differences phase. **\n")) ;
02526 
02527   /* for each column in pivot row */
02528   rp = &A [pivot_row_start] ;
02529   rp_end = rp + pivot_row_length ;
02530   while (rp < rp_end)
02531   {
02532       /* get a column */
02533       col = *rp++ ;
02534       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
02535       hash = 0 ;
02536       cur_score = 0 ;
02537       cp = &A [Col [col].start] ;
02538       /* compact the column */
02539       new_cp = cp ;
02540       cp_end = cp + Col [col].length ;
02541 
02542       DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
02543 
02544       while (cp < cp_end)
02545       {
02546     /* get a row */
02547     row = *cp++ ;
02548     ASSERT(row >= 0 && row < n_row) ;
02549     row_mark = Row [row].shared2.mark ;
02550     /* skip if dead */
02551     if (ROW_IS_MARKED_DEAD (row_mark))
02552     {
02553         DEBUG4 ((" Row %d, dead\n", row)) ;
02554         continue ;
02555     }
02556     DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
02557     ASSERT (row_mark >= tag_mark) ;
02558     /* compact the column */
02559     *new_cp++ = row ;
02560     /* compute hash function */
02561     hash += row ;
02562     /* add set difference */
02563     cur_score += row_mark - tag_mark ;
02564     /* integer overflow... */
02565     cur_score = MIN (cur_score, n_col) ;
02566       }
02567 
02568       /* recompute the column's length */
02569       Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
02570 
02571       /* === Further mass elimination ================================= */
02572 
02573       if (Col [col].length == 0)
02574       {
02575     DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
02576     /* nothing left but the pivot row in this column */
02577     KILL_PRINCIPAL_COL (col) ;
02578     pivot_row_degree -= Col [col].shared1.thickness ;
02579     ASSERT (pivot_row_degree >= 0) ;
02580     /* order it */
02581     Col [col].shared2.order = k ;
02582     /* increment order count by column thickness */
02583     k += Col [col].shared1.thickness ;
02584       }
02585       else
02586       {
02587     /* === Prepare for supercolumn detection ==================== */
02588 
02589     DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
02590 
02591     /* save score so far */
02592     Col [col].shared2.score = cur_score ;
02593 
02594     /* add column to hash table, for supercolumn detection */
02595     hash %= n_col + 1 ;
02596 
02597     DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
02598     ASSERT (((Int) hash) <= n_col) ;
02599 
02600     head_column = head [hash] ;
02601     if (head_column > EMPTY)
02602     {
02603         /* degree list "hash" is non-empty, use prev (shared3) of */
02604         /* first column in degree list as head of hash bucket */
02605         first_col = Col [head_column].shared3.headhash ;
02606         Col [head_column].shared3.headhash = col ;
02607     }
02608     else
02609     {
02610         /* degree list "hash" is empty, use head as hash bucket */
02611         first_col = - (head_column + 2) ;
02612         head [hash] = - (col + 2) ;
02613     }
02614     Col [col].shared4.hash_next = first_col ;
02615 
02616     /* save hash function in Col [col].shared3.hash */
02617     Col [col].shared3.hash = (Int) hash ;
02618     ASSERT (COL_IS_ALIVE (col)) ;
02619       }
02620   }
02621 
02622   /* The approximate external column degree is now computed.  */
02623 
02624   /* === Supercolumn detection ======================================== */
02625 
02626   DEBUG3 (("** Supercolumn detection phase. **\n")) ;
02627 
02628   detect_super_cols (
02629 
02630 #ifndef NDEBUG
02631     n_col, Row,
02632 #endif /* NDEBUG */
02633 
02634     Col, A, head, pivot_row_start, pivot_row_length) ;
02635 
02636   /* === Kill the pivotal column ====================================== */
02637 
02638   KILL_PRINCIPAL_COL (pivot_col) ;
02639 
02640   /* === Clear mark =================================================== */
02641 
02642   tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
02643 
02644 #ifndef NDEBUG
02645   DEBUG3 (("check3\n")) ;
02646   debug_mark (n_row, Row, tag_mark, max_mark) ;
02647 #endif /* NDEBUG */
02648 
02649   /* === Finalize the new pivot row, and column scores ================ */
02650 
02651   DEBUG3 (("** Finalize scores phase. **\n")) ;
02652 
02653   /* for each column in pivot row */
02654   rp = &A [pivot_row_start] ;
02655   /* compact the pivot row */
02656   new_rp = rp ;
02657   rp_end = rp + pivot_row_length ;
02658   while (rp < rp_end)
02659   {
02660       col = *rp++ ;
02661       /* skip dead columns */
02662       if (COL_IS_DEAD (col))
02663       {
02664     continue ;
02665       }
02666       *new_rp++ = col ;
02667       /* add new pivot row to column */
02668       A [Col [col].start + (Col [col].length++)] = pivot_row ;
02669 
02670       /* retrieve score so far and add on pivot row's degree. */
02671       /* (we wait until here for this in case the pivot */
02672       /* row's degree was reduced due to mass elimination). */
02673       cur_score = Col [col].shared2.score + pivot_row_degree ;
02674 
02675       /* calculate the max possible score as the number of */
02676       /* external columns minus the 'k' value minus the */
02677       /* columns thickness */
02678       max_score = n_col - k - Col [col].shared1.thickness ;
02679 
02680       /* make the score the external degree of the union-of-rows */
02681       cur_score -= Col [col].shared1.thickness ;
02682 
02683       /* make sure score is less or equal than the max score */
02684       cur_score = MIN (cur_score, max_score) ;
02685       ASSERT (cur_score >= 0) ;
02686 
02687       /* store updated score */
02688       Col [col].shared2.score = cur_score ;
02689 
02690       /* === Place column back in degree list ========================= */
02691 
02692       ASSERT (min_score >= 0) ;
02693       ASSERT (min_score <= n_col) ;
02694       ASSERT (cur_score >= 0) ;
02695       ASSERT (cur_score <= n_col) ;
02696       ASSERT (head [cur_score] >= EMPTY) ;
02697       next_col = head [cur_score] ;
02698       Col [col].shared4.degree_next = next_col ;
02699       Col [col].shared3.prev = EMPTY ;
02700       if (next_col != EMPTY)
02701       {
02702     Col [next_col].shared3.prev = col ;
02703       }
02704       head [cur_score] = col ;
02705 
02706       /* see if this score is less than current min */
02707       min_score = MIN (min_score, cur_score) ;
02708 
02709   }
02710 
02711 #ifndef NDEBUG
02712   debug_deg_lists (n_row, n_col, Row, Col, head,
02713     min_score, n_col2-k, max_deg) ;
02714 #endif /* NDEBUG */
02715 
02716   /* === Resurrect the new pivot row ================================== */
02717 
02718   if (pivot_row_degree > 0)
02719   {
02720       /* update pivot row length to reflect any cols that were killed */
02721       /* during super-col detection and mass elimination */
02722       Row [pivot_row].start  = pivot_row_start ;
02723       Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
02724       ASSERT (Row [pivot_row].length > 0) ;
02725       Row [pivot_row].shared1.degree = pivot_row_degree ;
02726       Row [pivot_row].shared2.mark = 0 ;
02727       /* pivot row is no longer dead */
02728 
02729       DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
02730       pivot_row, pivot_row_degree)) ;
02731   }
02732     }
02733 
02734     /* === All principal columns have now been ordered ====================== */
02735 
02736     return (ngarbage) ;
02737 }
02738 
02739 
02740 /* ========================================================================== */
02741 /* === order_children ======================================================= */
02742 /* ========================================================================== */
02743 
02744 /*
02745     The find_ordering routine has ordered all of the principal columns (the
02746     representatives of the supercolumns).  The non-principal columns have not
02747     yet been ordered.  This routine orders those columns by walking up the
02748     parent tree (a column is a child of the column which absorbed it).  The
02749     final permutation vector is then placed in p [0 ... n_col-1], with p [0]
02750     being the first column, and p [n_col-1] being the last.  It doesn't look
02751     like it at first glance, but be assured that this routine takes time linear
02752     in the number of columns.  Although not immediately obvious, the time
02753     taken by this routine is O (n_col), that is, linear in the number of
02754     columns.  Not user-callable.
02755 */
02756 
02757 PRIVATE void order_children
02758 (
02759     /* === Parameters ======================================================= */
02760 
02761     Int n_col,      /* number of columns of A */
02762     Colamd_Col Col [],    /* of size n_col+1 */
02763     Int p []      /* p [0 ... n_col-1] is the column permutation*/
02764 )
02765 {
02766     /* === Local variables ================================================== */
02767 
02768     Int i ;     /* loop counter for all columns */
02769     Int c ;     /* column index */
02770     Int parent ;    /* index of column's parent */
02771     Int order ;     /* column's order */
02772 
02773     /* === Order each non-principal column ================================== */
02774 
02775     for (i = 0 ; i < n_col ; i++)
02776     {
02777   /* find an un-ordered non-principal column */
02778   ASSERT (COL_IS_DEAD (i)) ;
02779   if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
02780   {
02781       parent = i ;
02782       /* once found, find its principal parent */
02783       do
02784       {
02785     parent = Col [parent].shared1.parent ;
02786       } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
02787 
02788       /* now, order all un-ordered non-principal columns along path */
02789       /* to this parent.  collapse tree at the same time */
02790       c = i ;
02791       /* get order of parent */
02792       order = Col [parent].shared2.order ;
02793 
02794       do
02795       {
02796     ASSERT (Col [c].shared2.order == EMPTY) ;
02797 
02798     /* order this column */
02799     Col [c].shared2.order = order++ ;
02800     /* collaps tree */
02801     Col [c].shared1.parent = parent ;
02802 
02803     /* get immediate parent of this column */
02804     c = Col [c].shared1.parent ;
02805 
02806     /* continue until we hit an ordered column.  There are */
02807     /* guarranteed not to be anymore unordered columns */
02808     /* above an ordered column */
02809       } while (Col [c].shared2.order == EMPTY) ;
02810 
02811       /* re-order the super_col parent to largest order for this group */
02812       Col [parent].shared2.order = order ;
02813   }
02814     }
02815 
02816     /* === Generate the permutation ========================================= */
02817 
02818     for (c = 0 ; c < n_col ; c++)
02819     {
02820   p [Col [c].shared2.order] = c ;
02821     }
02822 }
02823 
02824 
02825 /* ========================================================================== */
02826 /* === detect_super_cols ==================================================== */
02827 /* ========================================================================== */
02828 
02829 /*
02830     Detects supercolumns by finding matches between columns in the hash buckets.
02831     Check amongst columns in the set A [row_start ... row_start + row_length-1].
02832     The columns under consideration are currently *not* in the degree lists,
02833     and have already been placed in the hash buckets.
02834 
02835     The hash bucket for columns whose hash function is equal to h is stored
02836     as follows:
02837 
02838   if head [h] is >= 0, then head [h] contains a degree list, so:
02839 
02840     head [h] is the first column in degree bucket h.
02841     Col [head [h]].headhash gives the first column in hash bucket h.
02842 
02843   otherwise, the degree list is empty, and:
02844 
02845     -(head [h] + 2) is the first column in hash bucket h.
02846 
02847     For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
02848     column" pointer.  Col [c].shared3.hash is used instead as the hash number
02849     for that column.  The value of Col [c].shared4.hash_next is the next column
02850     in the same hash bucket.
02851 
02852     Assuming no, or "few" hash collisions, the time taken by this routine is
02853     linear in the sum of the sizes (lengths) of each column whose score has
02854     just been computed in the approximate degree computation.
02855     Not user-callable.
02856 */
02857 
02858 PRIVATE void detect_super_cols
02859 (
02860     /* === Parameters ======================================================= */
02861 
02862 #ifndef NDEBUG
02863     /* these two parameters are only needed when debugging is enabled: */
02864     Int n_col,      /* number of columns of A */
02865     Colamd_Row Row [],    /* of size n_row+1 */
02866 #endif /* NDEBUG */
02867 
02868     Colamd_Col Col [],    /* of size n_col+1 */
02869     Int A [],     /* row indices of A */
02870     Int head [],    /* head of degree lists and hash buckets */
02871     Int row_start,    /* pointer to set of columns to check */
02872     Int row_length    /* number of columns to check */
02873 )
02874 {
02875     /* === Local variables ================================================== */
02876 
02877     Int hash ;      /* hash value for a column */
02878     Int *rp ;     /* pointer to a row */
02879     Int c ;     /* a column index */
02880     Int super_c ;   /* column index of the column to absorb into */
02881     Int *cp1 ;      /* column pointer for column super_c */
02882     Int *cp2 ;      /* column pointer for column c */
02883     Int length ;    /* length of column super_c */
02884     Int prev_c ;    /* column preceding c in hash bucket */
02885     Int i ;     /* loop counter */
02886     Int *rp_end ;   /* pointer to the end of the row */
02887     Int col ;     /* a column index in the row to check */
02888     Int head_column ;   /* first column in hash bucket or degree list */
02889     Int first_col ;   /* first column in hash bucket */
02890 
02891     /* === Consider each column in the row ================================== */
02892 
02893     rp = &A [row_start] ;
02894     rp_end = rp + row_length ;
02895     while (rp < rp_end)
02896     {
02897   col = *rp++ ;
02898   if (COL_IS_DEAD (col))
02899   {
02900       continue ;
02901   }
02902 
02903   /* get hash number for this column */
02904   hash = Col [col].shared3.hash ;
02905   ASSERT (hash <= n_col) ;
02906 
02907   /* === Get the first column in this hash bucket ===================== */
02908 
02909   head_column = head [hash] ;
02910   if (head_column > EMPTY)
02911   {
02912       first_col = Col [head_column].shared3.headhash ;
02913   }
02914   else
02915   {
02916       first_col = - (head_column + 2) ;
02917   }
02918 
02919   /* === Consider each column in the hash bucket ====================== */
02920 
02921   for (super_c = first_col ; super_c != EMPTY ;
02922       super_c = Col [super_c].shared4.hash_next)
02923   {
02924       ASSERT (COL_IS_ALIVE (super_c)) ;
02925       ASSERT (Col [super_c].shared3.hash == hash) ;
02926       length = Col [super_c].length ;
02927 
02928       /* prev_c is the column preceding column c in the hash bucket */
02929       prev_c = super_c ;
02930 
02931       /* === Compare super_c with all columns after it ================ */
02932 
02933       for (c = Col [super_c].shared4.hash_next ;
02934      c != EMPTY ; c = Col [c].shared4.hash_next)
02935       {
02936     ASSERT (c != super_c) ;
02937     ASSERT (COL_IS_ALIVE (c)) ;
02938     ASSERT (Col [c].shared3.hash == hash) ;
02939 
02940     /* not identical if lengths or scores are different */
02941     if (Col [c].length != length ||
02942         Col [c].shared2.score != Col [super_c].shared2.score)
02943     {
02944         prev_c = c ;
02945         continue ;
02946     }
02947 
02948     /* compare the two columns */
02949     cp1 = &A [Col [super_c].start] ;
02950     cp2 = &A [Col [c].start] ;
02951 
02952     for (i = 0 ; i < length ; i++)
02953     {
02954         /* the columns are "clean" (no dead rows) */
02955         ASSERT (ROW_IS_ALIVE (*cp1))  ;
02956         ASSERT (ROW_IS_ALIVE (*cp2))  ;
02957         /* row indices will same order for both supercols, */
02958         /* no gather scatter nessasary */
02959         if (*cp1++ != *cp2++)
02960         {
02961       break ;
02962         }
02963     }
02964 
02965     /* the two columns are different if the for-loop "broke" */
02966     if (i != length)
02967     {
02968         prev_c = c ;
02969         continue ;
02970     }
02971 
02972     /* === Got it!  two columns are identical =================== */
02973 
02974     ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
02975 
02976     Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
02977     Col [c].shared1.parent = super_c ;
02978     KILL_NON_PRINCIPAL_COL (c) ;
02979     /* order c later, in order_children() */
02980     Col [c].shared2.order = EMPTY ;
02981     /* remove c from hash bucket */
02982     Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
02983       }
02984   }
02985 
02986   /* === Empty this hash bucket ======================================= */
02987 
02988   if (head_column > EMPTY)
02989   {
02990       /* corresponding degree list "hash" is not empty */
02991       Col [head_column].shared3.headhash = EMPTY ;
02992   }
02993   else
02994   {
02995       /* corresponding degree list "hash" is empty */
02996       head [hash] = EMPTY ;
02997   }
02998     }
02999 }
03000 
03001 
03002 /* ========================================================================== */
03003 /* === garbage_collection =================================================== */
03004 /* ========================================================================== */
03005 
03006 /*
03007     Defragments and compacts columns and rows in the workspace A.  Used when
03008     all avaliable memory has been used while performing row merging.  Returns
03009     the index of the first free position in A, after garbage collection.  The
03010     time taken by this routine is linear is the size of the array A, which is
03011     itself linear in the number of nonzeros in the input matrix.
03012     Not user-callable.
03013 */
03014 
03015 PRIVATE Int garbage_collection  /* returns the new value of pfree */
03016 (
03017     /* === Parameters ======================================================= */
03018 
03019     Int n_row,      /* number of rows */
03020     Int n_col,      /* number of columns */
03021     Colamd_Row Row [],    /* row info */
03022     Colamd_Col Col [],    /* column info */
03023     Int A [],     /* A [0 ... Alen-1] holds the matrix */
03024     Int *pfree      /* &A [0] ... pfree is in use */
03025 )
03026 {
03027     /* === Local variables ================================================== */
03028 
03029     Int *psrc ;     /* source pointer */
03030     Int *pdest ;    /* destination pointer */
03031     Int j ;     /* counter */
03032     Int r ;     /* a row index */
03033     Int c ;     /* a column index */
03034     Int length ;    /* length of a row or column */
03035 
03036 #ifndef NDEBUG
03037     Int debug_rows ;
03038     DEBUG2 (("Defrag..\n")) ;
03039     for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
03040     debug_rows = 0 ;
03041 #endif /* NDEBUG */
03042 
03043     /* === Defragment the columns =========================================== */
03044 
03045     pdest = &A[0] ;
03046     for (c = 0 ; c < n_col ; c++)
03047     {
03048   if (COL_IS_ALIVE (c))
03049   {
03050       psrc = &A [Col [c].start] ;
03051 
03052       /* move and compact the column */
03053       ASSERT (pdest <= psrc) ;
03054       Col [c].start = (Int) (pdest - &A [0]) ;
03055       length = Col [c].length ;
03056       for (j = 0 ; j < length ; j++)
03057       {
03058     r = *psrc++ ;
03059     if (ROW_IS_ALIVE (r))
03060     {
03061         *pdest++ = r ;
03062     }
03063       }
03064       Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
03065   }
03066     }
03067 
03068     /* === Prepare to defragment the rows =================================== */
03069 
03070     for (r = 0 ; r < n_row ; r++)
03071     {
03072   if (ROW_IS_DEAD (r) || (Row [r].length == 0))
03073   {
03074       /* This row is already dead, or is of zero length.  Cannot compact
03075        * a row of zero length, so kill it.  NOTE: in the current version,
03076        * there are no zero-length live rows.  Kill the row (for the first
03077        * time, or again) just to be safe. */
03078       KILL_ROW (r) ;
03079   }
03080   else
03081   {
03082       /* save first column index in Row [r].shared2.first_column */
03083       psrc = &A [Row [r].start] ;
03084       Row [r].shared2.first_column = *psrc ;
03085       ASSERT (ROW_IS_ALIVE (r)) ;
03086       /* flag the start of the row with the one's complement of row */
03087       *psrc = ONES_COMPLEMENT (r) ;
03088 #ifndef NDEBUG
03089       debug_rows++ ;
03090 #endif /* NDEBUG */
03091   }
03092     }
03093 
03094     /* === Defragment the rows ============================================== */
03095 
03096     psrc = pdest ;
03097     while (psrc < pfree)
03098     {
03099   /* find a negative number ... the start of a row */
03100   if (*psrc++ < 0)
03101   {
03102       psrc-- ;
03103       /* get the row index */
03104       r = ONES_COMPLEMENT (*psrc) ;
03105       ASSERT (r >= 0 && r < n_row) ;
03106       /* restore first column index */
03107       *psrc = Row [r].shared2.first_column ;
03108       ASSERT (ROW_IS_ALIVE (r)) ;
03109       ASSERT (Row [r].length > 0) ;
03110       /* move and compact the row */
03111       ASSERT (pdest <= psrc) ;
03112       Row [r].start = (Int) (pdest - &A [0]) ;
03113       length = Row [r].length ;
03114       for (j = 0 ; j < length ; j++)
03115       {
03116     c = *psrc++ ;
03117     if (COL_IS_ALIVE (c))
03118     {
03119         *pdest++ = c ;
03120     }
03121       }
03122       Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
03123       ASSERT (Row [r].length > 0) ;
03124 #ifndef NDEBUG
03125       debug_rows-- ;
03126 #endif /* NDEBUG */
03127   }
03128     }
03129     /* ensure we found all the rows */
03130     ASSERT (debug_rows == 0) ;
03131 
03132     /* === Return the new value of pfree ==================================== */
03133 
03134     return ((Int) (pdest - &A [0])) ;
03135 }
03136 
03137 
03138 /* ========================================================================== */
03139 /* === clear_mark =========================================================== */
03140 /* ========================================================================== */
03141 
03142 /*
03143     Clears the Row [].shared2.mark array, and returns the new tag_mark.
03144     Return value is the new tag_mark.  Not user-callable.
03145 */
03146 
03147 PRIVATE Int clear_mark  /* return the new value for tag_mark */
03148 (
03149     /* === Parameters ======================================================= */
03150 
03151     Int tag_mark, /* new value of tag_mark */
03152     Int max_mark, /* max allowed value of tag_mark */
03153 
03154     Int n_row,    /* number of rows in A */
03155     Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */
03156 )
03157 {
03158     /* === Local variables ================================================== */
03159 
03160     Int r ;
03161 
03162     if (tag_mark <= 0 || tag_mark >= max_mark)
03163     {
03164   for (r = 0 ; r < n_row ; r++)
03165   {
03166       if (ROW_IS_ALIVE (r))
03167       {
03168     Row [r].shared2.mark = 0 ;
03169       }
03170   }
03171   tag_mark = 1 ;
03172     }
03173 
03174     return (tag_mark) ;
03175 }
03176 
03177 
03178 /* ========================================================================== */
03179 /* === print_report ========================================================= */
03180 /* ========================================================================== */
03181 
03182 PRIVATE void print_report
03183 (
03184     char *method,
03185     Int stats [COLAMD_STATS]
03186 )
03187 {
03188 
03189     Int i1, i2, i3 ;
03190 
03191     PRINTF (("\n%s version %d.%d, %s: ", method,
03192       COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
03193 
03194     if (!stats)
03195     {
03196       PRINTF (("No statistics available.\n")) ;
03197   return ;
03198     }
03199 
03200     i1 = stats [COLAMD_INFO1] ;
03201     i2 = stats [COLAMD_INFO2] ;
03202     i3 = stats [COLAMD_INFO3] ;
03203 
03204     if (stats [COLAMD_STATUS] >= 0)
03205     {
03206       PRINTF (("OK.  ")) ;
03207     }
03208     else
03209     {
03210       PRINTF (("ERROR.  ")) ;
03211     }
03212 
03213     switch (stats [COLAMD_STATUS])
03214     {
03215 
03216   case COLAMD_OK_BUT_JUMBLED:
03217 
03218       PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
03219 
03220       PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
03221       method, i3)) ;
03222 
03223       PRINTF(("%s: last seen duplicate or out-of-order row index:   %d\n",
03224       method, INDEX (i2))) ;
03225 
03226       PRINTF(("%s: last seen in column:                             %d",
03227       method, INDEX (i1))) ;
03228 
03229       /* no break - fall through to next case instead */
03230 
03231   case COLAMD_OK:
03232 
03233       PRINTF(("\n")) ;
03234 
03235       PRINTF(("%s: number of dense or empty rows ignored:           %d\n",
03236       method, stats [COLAMD_DENSE_ROW])) ;
03237 
03238       PRINTF(("%s: number of dense or empty columns ignored:        %d\n",
03239       method, stats [COLAMD_DENSE_COL])) ;
03240 
03241       PRINTF(("%s: number of garbage collections performed:         %d\n",
03242       method, stats [COLAMD_DEFRAG_COUNT])) ;
03243       break ;
03244 
03245   case COLAMD_ERROR_A_not_present:
03246 
03247       PRINTF(("Array A (row indices of matrix) not present.\n")) ;
03248       break ;
03249 
03250   case COLAMD_ERROR_p_not_present:
03251 
03252       PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
03253       break ;
03254 
03255   case COLAMD_ERROR_nrow_negative:
03256 
03257       PRINTF(("Invalid number of rows (%d).\n", i1)) ;
03258       break ;
03259 
03260   case COLAMD_ERROR_ncol_negative:
03261 
03262       PRINTF(("Invalid number of columns (%d).\n", i1)) ;
03263       break ;
03264 
03265   case COLAMD_ERROR_nnz_negative:
03266 
03267       PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
03268       break ;
03269 
03270   case COLAMD_ERROR_p0_nonzero:
03271 
03272       PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
03273       break ;
03274 
03275   case COLAMD_ERROR_A_too_small:
03276 
03277       PRINTF(("Array A too small.\n")) ;
03278       PRINTF(("        Need Alen >= %d, but given only Alen = %d.\n",
03279       i1, i2)) ;
03280       break ;
03281 
03282   case COLAMD_ERROR_col_length_negative:
03283 
03284       PRINTF
03285       (("Column %d has a negative number of nonzero entries (%d).\n",
03286       INDEX (i1), i2)) ;
03287       break ;
03288 
03289   case COLAMD_ERROR_row_index_out_of_bounds:
03290 
03291       PRINTF
03292       (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
03293       INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
03294       break ;
03295 
03296   case COLAMD_ERROR_out_of_memory:
03297 
03298       PRINTF(("Out of memory.\n")) ;
03299       break ;
03300 
03301   /* v2.4: internal-error case deleted */
03302     }
03303 }
03304 
03305 
03306 
03307 
03308 /* ========================================================================== */
03309 /* === colamd debugging routines ============================================ */
03310 /* ========================================================================== */
03311 
03312 /* When debugging is disabled, the remainder of this file is ignored. */
03313 
03314 #ifndef NDEBUG
03315 
03316 
03317 /* ========================================================================== */
03318 /* === debug_structures ===================================================== */
03319 /* ========================================================================== */
03320 
03321 /*
03322     At this point, all empty rows and columns are dead.  All live columns
03323     are "clean" (containing no dead rows) and simplicial (no supercolumns
03324     yet).  Rows may contain dead columns, but all live rows contain at
03325     least one live column.
03326 */
03327 
03328 PRIVATE void debug_structures
03329 (
03330     /* === Parameters ======================================================= */
03331 
03332     Int n_row,
03333     Int n_col,
03334     Colamd_Row Row [],
03335     Colamd_Col Col [],
03336     Int A [],
03337     Int n_col2
03338 )
03339 {
03340     /* === Local variables ================================================== */
03341 
03342     Int i ;
03343     Int c ;
03344     Int *cp ;
03345     Int *cp_end ;
03346     Int len ;
03347     Int score ;
03348     Int r ;
03349     Int *rp ;
03350     Int *rp_end ;
03351     Int deg ;
03352 
03353     /* === Check A, Row, and Col ============================================ */
03354 
03355     for (c = 0 ; c < n_col ; c++)
03356     {
03357   if (COL_IS_ALIVE (c))
03358   {
03359       len = Col [c].length ;
03360       score = Col [c].shared2.score ;
03361       DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
03362       ASSERT (len > 0) ;
03363       ASSERT (score >= 0) ;
03364       ASSERT (Col [c].shared1.thickness == 1) ;
03365       cp = &A [Col [c].start] ;
03366       cp_end = cp + len ;
03367       while (cp < cp_end)
03368       {
03369     r = *cp++ ;
03370     ASSERT (ROW_IS_ALIVE (r)) ;
03371       }
03372   }
03373   else
03374   {
03375       i = Col [c].shared2.order ;
03376       ASSERT (i >= n_col2 && i < n_col) ;
03377   }
03378     }
03379 
03380     for (r = 0 ; r < n_row ; r++)
03381     {
03382   if (ROW_IS_ALIVE (r))
03383   {
03384       i = 0 ;
03385       len = Row [r].length ;
03386       deg = Row [r].shared1.degree ;
03387       ASSERT (len > 0) ;
03388       ASSERT (deg > 0) ;
03389       rp = &A [Row [r].start] ;
03390       rp_end = rp + len ;
03391       while (rp < rp_end)
03392       {
03393     c = *rp++ ;
03394     if (COL_IS_ALIVE (c))
03395     {
03396         i++ ;
03397     }
03398       }
03399       ASSERT (i > 0) ;
03400   }
03401     }
03402 }
03403 
03404 
03405 /* ========================================================================== */
03406 /* === debug_deg_lists ====================================================== */
03407 /* ========================================================================== */
03408 
03409 /*
03410     Prints the contents of the degree lists.  Counts the number of columns
03411     in the degree list and compares it to the total it should have.  Also
03412     checks the row degrees.
03413 */
03414 
03415 PRIVATE void debug_deg_lists
03416 (
03417     /* === Parameters ======================================================= */
03418 
03419     Int n_row,
03420     Int n_col,
03421     Colamd_Row Row [],
03422     Colamd_Col Col [],
03423     Int head [],
03424     Int min_score,
03425     Int should,
03426     Int max_deg
03427 )
03428 {
03429     /* === Local variables ================================================== */
03430 
03431     Int deg ;
03432     Int col ;
03433     Int have ;
03434     Int row ;
03435 
03436     /* === Check the degree lists =========================================== */
03437 
03438     if (n_col > 10000 && colamd_debug <= 0)
03439     {
03440   return ;
03441     }
03442     have = 0 ;
03443     DEBUG4 (("Degree lists: %d\n", min_score)) ;
03444     for (deg = 0 ; deg <= n_col ; deg++)
03445     {
03446   col = head [deg] ;
03447   if (col == EMPTY)
03448   {
03449       continue ;
03450   }
03451   DEBUG4 (("%d:", deg)) ;
03452   while (col != EMPTY)
03453   {
03454       DEBUG4 ((" %d", col)) ;
03455       have += Col [col].shared1.thickness ;
03456       ASSERT (COL_IS_ALIVE (col)) ;
03457       col = Col [col].shared4.degree_next ;
03458   }
03459   DEBUG4 (("\n")) ;
03460     }
03461     DEBUG4 (("should %d have %d\n", should, have)) ;
03462     ASSERT (should == have) ;
03463 
03464     /* === Check the row degrees ============================================ */
03465 
03466     if (n_row > 10000 && colamd_debug <= 0)
03467     {
03468   return ;
03469     }
03470     for (row = 0 ; row < n_row ; row++)
03471     {
03472   if (ROW_IS_ALIVE (row))
03473   {
03474       ASSERT (Row [row].shared1.degree <= max_deg) ;
03475   }
03476     }
03477 }
03478 
03479 
03480 /* ========================================================================== */
03481 /* === debug_mark =========================================================== */
03482 /* ========================================================================== */
03483 
03484 /*
03485     Ensures that the tag_mark is less that the maximum and also ensures that
03486     each entry in the mark array is less than the tag mark.
03487 */
03488 
03489 PRIVATE void debug_mark
03490 (
03491     /* === Parameters ======================================================= */
03492 
03493     Int n_row,
03494     Colamd_Row Row [],
03495     Int tag_mark,
03496     Int max_mark
03497 )
03498 {
03499     /* === Local variables ================================================== */
03500 
03501     Int r ;
03502 
03503     /* === Check the Row marks ============================================== */
03504 
03505     ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
03506     if (n_row > 10000 && colamd_debug <= 0)
03507     {
03508   return ;
03509     }
03510     for (r = 0 ; r < n_row ; r++)
03511     {
03512   ASSERT (Row [r].shared2.mark < tag_mark) ;
03513     }
03514 }
03515 
03516 
03517 /* ========================================================================== */
03518 /* === debug_matrix ========================================================= */
03519 /* ========================================================================== */
03520 
03521 /*
03522     Prints out the contents of the columns and the rows.
03523 */
03524 
03525 PRIVATE void debug_matrix
03526 (
03527     /* === Parameters ======================================================= */
03528 
03529     Int n_row,
03530     Int n_col,
03531     Colamd_Row Row [],
03532     Colamd_Col Col [],
03533     Int A []
03534 )
03535 {
03536     /* === Local variables ================================================== */
03537 
03538     Int r ;
03539     Int c ;
03540     Int *rp ;
03541     Int *rp_end ;
03542     Int *cp ;
03543     Int *cp_end ;
03544 
03545     /* === Dump the rows and columns of the matrix ========================== */
03546 
03547     if (colamd_debug < 3)
03548     {
03549   return ;
03550     }
03551     DEBUG3 (("DUMP MATRIX:\n")) ;
03552     for (r = 0 ; r < n_row ; r++)
03553     {
03554   DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ;
03555   if (ROW_IS_DEAD (r))
03556   {
03557       continue ;
03558   }
03559   DEBUG3 (("start %d length %d degree %d\n",
03560     Row [r].start, Row [r].length, Row [r].shared1.degree)) ;
03561   rp = &A [Row [r].start] ;
03562   rp_end = rp + Row [r].length ;
03563   while (rp < rp_end)
03564   {
03565       c = *rp++ ;
03566       DEBUG4 (("  %d col %d\n", COL_IS_ALIVE (c), c)) ;
03567   }
03568     }
03569 
03570     for (c = 0 ; c < n_col ; c++)
03571     {
03572   DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ;
03573   if (COL_IS_DEAD (c))
03574   {
03575       continue ;
03576   }
03577   DEBUG3 (("start %d length %d shared1 %d shared2 %d\n",
03578     Col [c].start, Col [c].length,
03579     Col [c].shared1.thickness, Col [c].shared2.score)) ;
03580   cp = &A [Col [c].start] ;
03581   cp_end = cp + Col [c].length ;
03582   while (cp < cp_end)
03583   {
03584       r = *cp++ ;
03585       DEBUG4 (("  %d row %d\n", ROW_IS_ALIVE (r), r)) ;
03586   }
03587     }
03588 }
03589 
03590 PRIVATE void colamd_get_debug
03591 (
03592     char *method
03593 )
03594 {
03595     FILE *f ;
03596     colamd_debug = 0 ;    /* no debug printing */
03597     f = fopen ("debug", "r") ;
03598     if (f == (FILE *) NULL)
03599     {
03600   colamd_debug = 0 ;
03601     }
03602     else
03603     {
03604   fscanf (f, "%d", &colamd_debug) ;
03605   fclose (f) ;
03606     }
03607     DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
03608       method, colamd_debug)) ;
03609 }
03610 
03611 #endif /* NDEBUG */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines