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