extern/colamd/Source/colamd.c

   1 /* ========================================================================== */
   2 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */
   3 /* ========================================================================== */
   4
   5 /* COLAMD / SYMAMD
   6
   7     colamd:  an approximate minimum degree column ordering algorithm,
   8         for LU factorization of symmetric or unsymmetric matrices,
   9         QR factorization, least squares, interior point methods for
  10         linear programming problems, and other related problems.
  11
  12     symamd:  an approximate minimum degree ordering algorithm for Cholesky
  13         factorization of symmetric matrices.
  14
  15     Purpose:
  16
  17         Colamd computes a permutation Q such that the Cholesky factorization of
  18         (AQ)'(AQ) has less fill-in and requires fewer floating point operations
  19         than A'A.  This also provides a good ordering for sparse partial
  20         pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
  21         factorization, and P is computed during numerical factorization via
  22         conventional partial pivoting with row interchanges.  Colamd is the
  23         column ordering method used in SuperLU, part of the ScaLAPACK library.
  24         It is also available as built-in function in MATLAB Version 6,
  25         available from MathWorks, Inc. (http://www.mathworks.com).  This
  26         routine can be used in place of colmmd in MATLAB.
  27
  28         Symamd computes a permutation P of a symmetric matrix A such that the
  29         Cholesky factorization of PAP' has less fill-in and requires fewer
  30         floating point operations than A.  Symamd constructs a matrix M such
  31         that M'M has the same nonzero pattern of A, and then orders the columns
  32         of M using colmmd.  The column ordering of M is then returned as the
  33         row and column ordering P of A.
  34
  35     Authors:
  36
  37         The authors of the code itself are Stefan I. Larimore and Timothy A.
  38         Davis (davis at cise.ufl.edu), University of Florida.  The algorithm was
  39         developed in collaboration with John Gilbert, Xerox PARC, and Esmond
  40         Ng, Oak Ridge National Laboratory.
  41
  42     Acknowledgements:
  43
  44         This work was supported by the National Science Foundation, under
  45         grants DMS-9504974 and DMS-9803599.
  46
  47     Copyright and License:
  48
  49         Copyright (c) 1998-2007, Timothy A. Davis, All Rights Reserved.
  50         COLAMD is also available under alternate licenses, contact T. Davis
  51         for details.
  52
  53         This library is free software; you can redistribute it and/or
  54         modify it under the terms of the GNU Lesser General Public
  55         License as published by the Free Software Foundation; either
  56         version 2.1 of the License, or (at your option) any later version.
  57
  58         This library is distributed in the hope that it will be useful,
  59         but WITHOUT ANY WARRANTY; without even the implied warranty of
  60         MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  61         Lesser General Public License for more details.
  62
  63         You should have received a copy of the GNU Lesser General Public
  64         License along with this library; if not, write to the Free Software
  65         Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
  66         USA
  67
  68         Permission is hereby granted to use or copy this program under the
  69         terms of the GNU LGPL, provided that the Copyright, this License,
  70         and the Availability of the original version is retained on all copies.
  71         User documentation of any code that uses this code or any modified
  72         version of this code must cite the Copyright, this License, the
  73         Availability note, and "Used by permission." Permission to modify
  74         the code and to distribute modified code is granted, provided the
  75         Copyright, this License, and the Availability note are retained,
  76         and a notice that the code was modified is included.
  77
  78     Availability:
  79
  80         The colamd/symamd library is available at
  81
  82             http://www.cise.ufl.edu/research/sparse/colamd/
  83
  84         This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c
  85         file.  It requires the colamd.h file.  It is required by the colamdmex.c
  86         and symamdmex.c files, for the MATLAB interface to colamd and symamd.
  87         Appears as ACM Algorithm 836.
  88
  89     See the ChangeLog file for changes since Version 1.0.
  90
  91     References:
  92
  93         T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, An approximate column
  94         minimum degree ordering algorithm, ACM Transactions on Mathematical
  95         Software, vol. 30, no. 3., pp. 353-376, 2004.
  96
  97         T. A. Davis, J. R. Gilbert, S. Larimore, E. Ng, Algorithm 836: COLAMD,
  98         an approximate column minimum degree ordering algorithm, ACM
  99         Transactions on Mathematical Software, vol. 30, no. 3., pp. 377-380,
 100         2004.
 101
 102 */
 103
 104 /* ========================================================================== */
 105 /* === Description of user-callable routines ================================ */
 106 /* ========================================================================== */
 107
 108 /* COLAMD includes both int and UF_long versions of all its routines.  The
 109  * description below is for the int version.  For UF_long, all int arguments
 110  * become UF_long.  UF_long is normally defined as long, except for WIN64.
 111
 112     ----------------------------------------------------------------------------
 113     colamd_recommended:
 114     ----------------------------------------------------------------------------
 115
 116         C syntax:
 117
 118             #include "colamd.h"
 119             size_t colamd_recommended (int nnz, int n_row, int n_col) ;
 120             size_t colamd_l_recommended (UF_long nnz, UF_long n_row,
 121                 UF_long n_col) ;
 122
 123         Purpose:
 124
 125             Returns recommended value of Alen for use by colamd.  Returns 0
 126             if any input argument is negative.  The use of this routine
 127             is optional.  Not needed for symamd, which dynamically allocates
 128             its own memory.
 129
 130             Note that in v2.4 and earlier, these routines returned int or long.
 131             They now return a value of type size_t.
 132
 133         Arguments (all input arguments):
 134
 135             int nnz ;           Number of nonzeros in the matrix A.  This must
 136                                 be the same value as p [n_col] in the call to
 137                                 colamd - otherwise you will get a wrong value
 138                                 of the recommended memory to use.
 139
 140             int n_row ;         Number of rows in the matrix A.
 141
 142             int n_col ;         Number of columns in the matrix A.
 143
 144     ----------------------------------------------------------------------------
 145     colamd_set_defaults:
 146     ----------------------------------------------------------------------------
 147
 148         C syntax:
 149
 150             #include "colamd.h"
 151             colamd_set_defaults (double knobs [COLAMD_KNOBS]) ;
 152             colamd_l_set_defaults (double knobs [COLAMD_KNOBS]) ;
 153
 154         Purpose:
 155
 156             Sets the default parameters.  The use of this routine is optional.
 157
 158         Arguments:
 159
 160             double knobs [COLAMD_KNOBS] ;       Output only.
 161
 162                 NOTE: the meaning of the dense row/col knobs has changed in v2.4
 163
 164                 knobs [0] and knobs [1] control dense row and col detection:
 165
 166                 Colamd: rows with more than
 167                 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n_col))
 168                 entries are removed prior to ordering.  Columns with more than
 169                 max (16, knobs [COLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
 170                 entries are removed prior to
 171                 ordering, and placed last in the output column ordering.
 172
 173                 Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0].
 174                 Rows and columns with more than
 175                 max (16, knobs [COLAMD_DENSE_ROW] * sqrt (n))
 176                 entries are removed prior to ordering, and placed last in the
 177                 output ordering.
 178
 179                 COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1,
 180                 respectively, in colamd.h.  Default values of these two knobs
 181                 are both 10.  Currently, only knobs [0] and knobs [1] are
 182                 used, but future versions may use more knobs.  If so, they will
 183                 be properly set to their defaults by the future version of
 184                 colamd_set_defaults, so that the code that calls colamd will
 185                 not need to change, assuming that you either use
 186                 colamd_set_defaults, or pass a (double *) NULL pointer as the
 187                 knobs array to colamd or symamd.
 188
 189             knobs [2]: aggressive absorption
 190
 191                 knobs [COLAMD_AGGRESSIVE] controls whether or not to do
 192                 aggressive absorption during the ordering.  Default is TRUE.
 193
 194
 195     ----------------------------------------------------------------------------
 196     colamd:
 197     ----------------------------------------------------------------------------
 198
 199         C syntax:
 200
 201             #include "colamd.h"
 202             int colamd (int n_row, int n_col, int Alen, int *A, int *p,
 203                 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ;
 204             UF_long colamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
 205                 UF_long *A, UF_long *p, double knobs [COLAMD_KNOBS],
 206                 UF_long stats [COLAMD_STATS]) ;
 207
 208         Purpose:
 209
 210             Computes a column ordering (Q) of A such that P(AQ)=LU or
 211             (AQ)'AQ=LL' have less fill-in and require fewer floating point
 212             operations than factorizing the unpermuted matrix A or A'A,
 213             respectively.
 214
 215         Returns:
 216
 217             TRUE (1) if successful, FALSE (0) otherwise.
 218
 219         Arguments:
 220
 221             int n_row ;         Input argument.
 222
 223                 Number of rows in the matrix A.
 224                 Restriction:  n_row >= 0.
 225                 Colamd returns FALSE if n_row is negative.
 226
 227             int n_col ;         Input argument.
 228
 229                 Number of columns in the matrix A.
 230                 Restriction:  n_col >= 0.
 231                 Colamd returns FALSE if n_col is negative.
 232
 233             int Alen ;          Input argument.
 234
 235                 Restriction (see note):
 236                 Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col
 237                 Colamd returns FALSE if these conditions are not met.
 238
 239                 Note:  this restriction makes an modest assumption regarding
 240                 the size of the two typedef's structures in colamd.h.
 241                 We do, however, guarantee that
 242
 243                         Alen >= colamd_recommended (nnz, n_row, n_col)
 244
 245                 will be sufficient.  Note: the macro version does not check
 246                 for integer overflow, and thus is not recommended.  Use
 247                 the colamd_recommended routine instead.
 248
 249             int A [Alen] ;      Input argument, undefined on output.
 250
 251                 A is an integer array of size Alen.  Alen must be at least as
 252                 large as the bare minimum value given above, but this is very
 253                 low, and can result in excessive run time.  For best
 254                 performance, we recommend that Alen be greater than or equal to
 255                 colamd_recommended (nnz, n_row, n_col), which adds
 256                 nnz/5 to the bare minimum value given above.
 257
 258                 On input, the row indices of the entries in column c of the
 259                 matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
 260                 in a given column c need not be in ascending order, and
 261                 duplicate row indices may be be present.  However, colamd will
 262                 work a little faster if both of these conditions are met
 263                 (Colamd puts the matrix into this format, if it finds that the
 264                 the conditions are not met).
 265
 266                 The matrix is 0-based.  That is, rows are in the range 0 to
 267                 n_row-1, and columns are in the range 0 to n_col-1.  Colamd
 268                 returns FALSE if any row index is out of range.
 269
 270                 The contents of A are modified during ordering, and are
 271                 undefined on output.
 272
 273             int p [n_col+1] ;   Both input and output argument.
 274
 275                 p is an integer array of size n_col+1.  On input, it holds the
 276                 "pointers" for the column form of the matrix A.  Column c of
 277                 the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
 278                 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
 279                 for all c in the range 0 to n_col-1.  The value p [n_col] is
 280                 thus the total number of entries in the pattern of the matrix A.
 281                 Colamd returns FALSE if these conditions are not met.
 282
 283                 On output, if colamd returns TRUE, the array p holds the column
 284                 permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
 285                 the first column index in the new ordering, and p [n_col-1] is
 286                 the last.  That is, p [k] = j means that column j of A is the
 287                 kth pivot column, in AQ, where k is in the range 0 to n_col-1
 288                 (p [0] = j means that column j of A is the first column in AQ).
 289
 290                 If colamd returns FALSE, then no permutation is returned, and
 291                 p is undefined on output.
 292
 293             double knobs [COLAMD_KNOBS] ;       Input argument.
 294
 295                 See colamd_set_defaults for a description.
 296
 297             int stats [COLAMD_STATS] ;          Output argument.
 298
 299                 Statistics on the ordering, and error status.
 300                 See colamd.h for related definitions.
 301                 Colamd returns FALSE if stats is not present.
 302
 303                 stats [0]:  number of dense or empty rows ignored.
 304
 305                 stats [1]:  number of dense or empty columns ignored (and
 306                                 ordered last in the output permutation p)
 307                                 Note that a row can become "empty" if it
 308                                 contains only "dense" and/or "empty" columns,
 309                                 and similarly a column can become "empty" if it
 310                                 only contains "dense" and/or "empty" rows.
 311
 312                 stats [2]:  number of garbage collections performed.
 313                                 This can be excessively high if Alen is close
 314                                 to the minimum required value.
 315
 316                 stats [3]:  status code.  < 0 is an error code.
 317                             > 1 is a warning or notice.
 318
 319                         0       OK.  Each column of the input matrix contained
 320                                 row indices in increasing order, with no
 321                                 duplicates.
 322
 323                         1       OK, but columns of input matrix were jumbled
 324                                 (unsorted columns or duplicate entries).  Colamd
 325                                 had to do some extra work to sort the matrix
 326                                 first and remove duplicate entries, but it
 327                                 still was able to return a valid permutation
 328                                 (return value of colamd was TRUE).
 329
 330                                         stats [4]: highest numbered column that
 331                                                 is unsorted or has duplicate
 332                                                 entries.
 333                                         stats [5]: last seen duplicate or
 334                                                 unsorted row index.
 335                                         stats [6]: number of duplicate or
 336                                                 unsorted row indices.
 337
 338                         -1      A is a null pointer
 339
 340                         -2      p is a null pointer
 341
 342                         -3      n_row is negative
 343
 344                                         stats [4]: n_row
 345
 346                         -4      n_col is negative
 347
 348                                         stats [4]: n_col
 349
 350                         -5      number of nonzeros in matrix is negative
 351
 352                                         stats [4]: number of nonzeros, p [n_col]
 353
 354                         -6      p [0] is nonzero
 355
 356                                         stats [4]: p [0]
 357
 358                         -7      A is too small
 359
 360                                         stats [4]: required size
 361                                         stats [5]: actual size (Alen)
 362
 363                         -8      a column has a negative number of entries
 364
 365                                         stats [4]: column with < 0 entries
 366                                         stats [5]: number of entries in col
 367
 368                         -9      a row index is out of bounds
 369
 370                                         stats [4]: column with bad row index
 371                                         stats [5]: bad row index
 372                                         stats [6]: n_row, # of rows of matrx
 373
 374                         -10     (unused; see symamd.c)
 375
 376                         -999    (unused; see symamd.c)
 377
 378                 Future versions may return more statistics in the stats array.
 379
 380         Example:
 381
 382             See http://www.cise.ufl.edu/research/sparse/colamd/example.c
 383             for a complete example.
 384
 385             To order the columns of a 5-by-4 matrix with 11 nonzero entries in
 386             the following nonzero pattern
 387
 388                 x 0 x 0
 389                 x 0 x x
 390                 0 x x 0
 391                 0 0 x x
 392                 x x 0 0
 393
 394             with default knobs and no output statistics, do the following:
 395
 396                 #include "colamd.h"
 397                 #define ALEN 100
 398                 int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
 399                 int p [ ] = {0, 3, 5, 9, 11} ;
 400                 int stats [COLAMD_STATS] ;
 401                 colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ;
 402
 403             The permutation is returned in the array p, and A is destroyed.
 404
 405     ----------------------------------------------------------------------------
 406     symamd:
 407     ----------------------------------------------------------------------------
 408
 409         C syntax:
 410
 411             #include "colamd.h"
 412             int symamd (int n, int *A, int *p, int *perm,
 413                 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS],
 414                 void (*allocate) (size_t, size_t), void (*release) (void *)) ;
 415             UF_long symamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
 416                 double knobs [COLAMD_KNOBS], UF_long stats [COLAMD_STATS],
 417                 void (*allocate) (size_t, size_t), void (*release) (void *)) ;
 418
 419         Purpose:
 420
 421             The symamd routine computes an ordering P of a symmetric sparse
 422             matrix A such that the Cholesky factorization PAP' = LL' remains
 423             sparse.  It is based on a column ordering of a matrix M constructed
 424             so that the nonzero pattern of M'M is the same as A.  The matrix A
 425             is assumed to be symmetric; only the strictly lower triangular part
 426             is accessed.  You must pass your selected memory allocator (usually
 427             calloc/free or mxCalloc/mxFree) to symamd, for it to allocate
 428             memory for the temporary matrix M.
 429
 430         Returns:
 431
 432             TRUE (1) if successful, FALSE (0) otherwise.
 433
 434         Arguments:
 435
 436             int n ;             Input argument.
 437
 438                 Number of rows and columns in the symmetrix matrix A.
 439                 Restriction:  n >= 0.
 440                 Symamd returns FALSE if n is negative.
 441
 442             int A [nnz] ;       Input argument.
 443
 444                 A is an integer array of size nnz, where nnz = p [n].
 445
 446                 The row indices of the entries in column c of the matrix are
 447                 held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
 448                 given column c need not be in ascending order, and duplicate
 449                 row indices may be present.  However, symamd will run faster
 450                 if the columns are in sorted order with no duplicate entries.
 451
 452                 The matrix is 0-based.  That is, rows are in the range 0 to
 453                 n-1, and columns are in the range 0 to n-1.  Symamd
 454                 returns FALSE if any row index is out of range.
 455
 456                 The contents of A are not modified.
 457
 458             int p [n+1] ;       Input argument.
 459
 460                 p is an integer array of size n+1.  On input, it holds the
 461                 "pointers" for the column form of the matrix A.  Column c of
 462                 the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
 463                 entry, p [0], must be zero, and p [c] <= p [c+1] must hold
 464                 for all c in the range 0 to n-1.  The value p [n] is
 465                 thus the total number of entries in the pattern of the matrix A.
 466                 Symamd returns FALSE if these conditions are not met.
 467
 468                 The contents of p are not modified.
 469
 470             int perm [n+1] ;    Output argument.
 471
 472                 On output, if symamd returns TRUE, the array perm holds the
 473                 permutation P, where perm [0] is the first index in the new
 474                 ordering, and perm [n-1] is the last.  That is, perm [k] = j
 475                 means that row and column j of A is the kth column in PAP',
 476                 where k is in the range 0 to n-1 (perm [0] = j means
 477                 that row and column j of A are the first row and column in
 478                 PAP').  The array is used as a workspace during the ordering,
 479                 which is why it must be of length n+1, not just n.
 480
 481             double knobs [COLAMD_KNOBS] ;       Input argument.
 482
 483                 See colamd_set_defaults for a description.
 484
 485             int stats [COLAMD_STATS] ;          Output argument.
 486
 487                 Statistics on the ordering, and error status.
 488                 See colamd.h for related definitions.
 489                 Symamd returns FALSE if stats is not present.
 490
 491                 stats [0]:  number of dense or empty row and columns ignored
 492                                 (and ordered last in the output permutation
 493                                 perm).  Note that a row/column can become
 494                                 "empty" if it contains only "dense" and/or
 495                                 "empty" columns/rows.
 496
 497                 stats [1]:  (same as stats [0])
 498
 499                 stats [2]:  number of garbage collections performed.
 500
 501                 stats [3]:  status code.  < 0 is an error code.
 502                             > 1 is a warning or notice.
 503
 504                         0       OK.  Each column of the input matrix contained
 505                                 row indices in increasing order, with no
 506                                 duplicates.
 507
 508                         1       OK, but columns of input matrix were jumbled
 509                                 (unsorted columns or duplicate entries).  Symamd
 510                                 had to do some extra work to sort the matrix
 511                                 first and remove duplicate entries, but it
 512                                 still was able to return a valid permutation
 513                                 (return value of symamd was TRUE).
 514
 515                                         stats [4]: highest numbered column that
 516                                                 is unsorted or has duplicate
 517                                                 entries.
 518                                         stats [5]: last seen duplicate or
 519                                                 unsorted row index.
 520                                         stats [6]: number of duplicate or
 521                                                 unsorted row indices.
 522
 523                         -1      A is a null pointer
 524
 525                         -2      p is a null pointer
 526
 527                         -3      (unused, see colamd.c)
 528
 529                         -4      n is negative
 530
 531                                         stats [4]: n
 532
 533                         -5      number of nonzeros in matrix is negative
 534
 535                                         stats [4]: # of nonzeros (p [n]).
 536
 537                         -6      p [0] is nonzero
 538
 539                                         stats [4]: p [0]
 540
 541                         -7      (unused)
 542
 543                         -8      a column has a negative number of entries
 544
 545                                         stats [4]: column with < 0 entries
 546                                         stats [5]: number of entries in col
 547
 548                         -9      a row index is out of bounds
 549
 550                                         stats [4]: column with bad row index
 551                                         stats [5]: bad row index
 552                                         stats [6]: n_row, # of rows of matrx
 553
 554                         -10     out of memory (unable to allocate temporary
 555                                 workspace for M or count arrays using the
 556                                 "allocate" routine passed into symamd).
 557
 558                 Future versions may return more statistics in the stats array.
 559
 560             void * (*allocate) (size_t, size_t)
 561
 562                 A pointer to a function providing memory allocation.  The
 563                 allocated memory must be returned initialized to zero.  For a
 564                 C application, this argument should normally be a pointer to
 565                 calloc.  For a MATLAB mexFunction, the routine mxCalloc is
 566                 passed instead.
 567
 568             void (*release) (size_t, size_t)
 569
 570                 A pointer to a function that frees memory allocated by the
 571                 memory allocation routine above.  For a C application, this
 572                 argument should normally be a pointer to free.  For a MATLAB
 573                 mexFunction, the routine mxFree is passed instead.
 574
 575
 576     ----------------------------------------------------------------------------
 577     colamd_report:
 578     ----------------------------------------------------------------------------
 579
 580         C syntax:
 581
 582             #include "colamd.h"
 583             colamd_report (int stats [COLAMD_STATS]) ;
 584             colamd_l_report (UF_long stats [COLAMD_STATS]) ;
 585
 586         Purpose:
 587
 588             Prints the error status and statistics recorded in the stats
 589             array on the standard error output (for a standard C routine)
 590             or on the MATLAB output (for a mexFunction).
 591
 592         Arguments:
 593
 594             int stats [COLAMD_STATS] ;  Input only.  Statistics from colamd.
 595
 596
 597     ----------------------------------------------------------------------------
 598     symamd_report:
 599     ----------------------------------------------------------------------------
 600
 601         C syntax:
 602
 603             #include "colamd.h"
 604             symamd_report (int stats [COLAMD_STATS]) ;
 605             symamd_l_report (UF_long stats [COLAMD_STATS]) ;
 606
 607         Purpose:
 608
 609             Prints the error status and statistics recorded in the stats
 610             array on the standard error output (for a standard C routine)
 611             or on the MATLAB output (for a mexFunction).
 612
 613         Arguments:
 614
 615             int stats [COLAMD_STATS] ;  Input only.  Statistics from symamd.
 616
 617
 618 */
 619
 620 /* ========================================================================== */
 621 /* === Scaffolding code definitions  ======================================== */
 622 /* ========================================================================== */
 623
 624 /* Ensure that debugging is turned off: */
 625 #ifndef NDEBUG
 626 #define NDEBUG
 627 #endif
 628
 629 /* turn on debugging by uncommenting the following line
 630  #undef NDEBUG
 631 */
 632
 633 /*
 634    Our "scaffolding code" philosophy:  In our opinion, well-written library
 635    code should keep its "debugging" code, and just normally have it turned off
 636    by the compiler so as not to interfere with performance.  This serves
 637    several purposes:
 638
 639    (1) assertions act as comments to the reader, telling you what the code
 640         expects at that point.  All assertions will always be true (unless
 641         there really is a bug, of course).
 642
 643    (2) leaving in the scaffolding code assists anyone who would like to modify
 644         the code, or understand the algorithm (by reading the debugging output,
 645         one can get a glimpse into what the code is doing).
 646
 647    (3) (gasp!) for actually finding bugs.  This code has been heavily tested
 648         and "should" be fully functional and bug-free ... but you never know...
 649
 650     The code will become outrageously slow when debugging is
 651     enabled.  To control the level of debugging output, set an environment
 652     variable D to 0 (little), 1 (some), 2, 3, or 4 (lots).  When debugging,
 653     you should see the following message on the standard output:
 654
 655         colamd: debug version, D = 1 (THIS WILL BE SLOW!)
 656
 657     or a similar message for symamd.  If you don't, then debugging has not
 658     been enabled.
 659
 660 */
 661
 662 /* ========================================================================== */
 663 /* === Include files ======================================================== */
 664 /* ========================================================================== */
 665
 666 #include "colamd.h"
 667 #include <limits.h>
 668 #include <math.h>
 669
 670 #ifdef MATLAB_MEX_FILE
 671 #include "mex.h"
 672 #include "matrix.h"
 673 #endif /* MATLAB_MEX_FILE */
 674
 675 #if !defined (NPRINT) || !defined (NDEBUG)
 676 #include <stdio.h>
 677 #endif
 678
 679 #ifndef NULL
 680 #define NULL ((void *) 0)
 681 #endif
 682
 683 /* ========================================================================== */
 684 /* === int or UF_long ======================================================= */
 685 /* ========================================================================== */
 686
 687 /* define UF_long */
 688 #include "UFconfig.h"
 689
 690 #ifdef DLONG
 691
 692 #define Int UF_long
 693 #define ID  UF_long_id
 694 #define Int_MAX UF_long_max
 695
 696 #define COLAMD_recommended colamd_l_recommended
 697 #define COLAMD_set_defaults colamd_l_set_defaults
 698 #define COLAMD_MAIN colamd_l
 699 #define SYMAMD_MAIN symamd_l
 700 #define COLAMD_report colamd_l_report
 701 #define SYMAMD_report symamd_l_report
 702
 703 #else
 704
 705 #define Int int
 706 #define ID "%d"
 707 #define Int_MAX INT_MAX
 708
 709 #define COLAMD_recommended colamd_recommended
 710 #define COLAMD_set_defaults colamd_set_defaults
 711 #define COLAMD_MAIN colamd
 712 #define SYMAMD_MAIN symamd
 713 #define COLAMD_report colamd_report
 714 #define SYMAMD_report symamd_report
 715
 716 #endif
 717
 718 /* ========================================================================== */
 719 /* === Row and Column structures ============================================ */
 720 /* ========================================================================== */
 721
 722 /* User code that makes use of the colamd/symamd routines need not directly */
 723 /* reference these structures.  They are used only for colamd_recommended. */
 724
 725 typedef struct Colamd_Col_struct
 726 {
 727     Int start ;         /* index for A of first row in this column, or DEAD */
 728                         /* if column is dead */
 729     Int length ;        /* number of rows in this column */
 730     union
 731     {
 732         Int thickness ; /* number of original columns represented by this */
 733                         /* col, if the column is alive */
 734         Int parent ;    /* parent in parent tree super-column structure, if */
 735                         /* the column is dead */
 736     } shared1 ;
 737     union
 738     {
 739         Int score ;     /* the score used to maintain heap, if col is alive */
 740         Int order ;     /* pivot ordering of this column, if col is dead */
 741     } shared2 ;
 742     union
 743     {
 744         Int headhash ;  /* head of a hash bucket, if col is at the head of */
 745                         /* a degree list */
 746         Int hash ;      /* hash value, if col is not in a degree list */
 747         Int prev ;      /* previous column in degree list, if col is in a */
 748                         /* degree list (but not at the head of a degree list) */
 749     } shared3 ;
 750     union
 751     {
 752         Int degree_next ;       /* next column, if col is in a degree list */
 753         Int hash_next ;         /* next column, if col is in a hash list */
 754     } shared4 ;
 755
 756 } Colamd_Col ;
 757
 758 typedef struct Colamd_Row_struct
 759 {
 760     Int start ;         /* index for A of first col in this row */
 761     Int length ;        /* number of principal columns in this row */
 762     union
 763     {
 764         Int degree ;    /* number of principal & non-principal columns in row */
 765         Int p ;         /* used as a row pointer in init_rows_cols () */
 766     } shared1 ;
 767     union
 768     {
 769         Int mark ;      /* for computing set differences and marking dead rows*/
 770         Int first_column ;/* first column in row (used in garbage collection) */
 771     } shared2 ;
 772
 773 } Colamd_Row ;
 774
 775 /* ========================================================================== */
 776 /* === Definitions ========================================================== */
 777 /* ========================================================================== */
 778
 779 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
 780 #define PUBLIC
 781 #define PRIVATE static
 782
 783 #define DENSE_DEGREE(alpha,n) \
 784     ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
 785
 786 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
 787 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
 788
 789 #define ONES_COMPLEMENT(r) (-(r)-1)
 790
 791 /* -------------------------------------------------------------------------- */
 792 /* Change for version 2.1:  define TRUE and FALSE only if not yet defined */
 793 /* -------------------------------------------------------------------------- */
 794
 795 #ifndef TRUE
 796 #define TRUE (1)
 797 #endif
 798
 799 #ifndef FALSE
 800 #define FALSE (0)
 801 #endif
 802
 803 /* -------------------------------------------------------------------------- */
 804
 805 #define EMPTY   (-1)
 806
 807 /* Row and column status */
 808 #define ALIVE   (0)
 809 #define DEAD    (-1)
 810
 811 /* Column status */
 812 #define DEAD_PRINCIPAL          (-1)
 813 #define DEAD_NON_PRINCIPAL      (-2)
 814
 815 /* Macros for row and column status update and checking. */
 816 #define ROW_IS_DEAD(r)                  ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
 817 #define ROW_IS_MARKED_DEAD(row_mark)    (row_mark < ALIVE)
 818 #define ROW_IS_ALIVE(r)                 (Row [r].shared2.mark >= ALIVE)
 819 #define COL_IS_DEAD(c)                  (Col [c].start < ALIVE)
 820 #define COL_IS_ALIVE(c)                 (Col [c].start >= ALIVE)
 821 #define COL_IS_DEAD_PRINCIPAL(c)        (Col [c].start == DEAD_PRINCIPAL)
 822 #define KILL_ROW(r)                     { Row [r].shared2.mark = DEAD ; }
 823 #define KILL_PRINCIPAL_COL(c)           { Col [c].start = DEAD_PRINCIPAL ; }
 824 #define KILL_NON_PRINCIPAL_COL(c)       { Col [c].start = DEAD_NON_PRINCIPAL ; }
 825
 826 /* ========================================================================== */
 827 /* === Colamd reporting mechanism =========================================== */
 828 /* ========================================================================== */
 829
 830 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
 831 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
 832 #define INDEX(i) ((i)+1)
 833 #else
 834 /* In C, matrices are 0-based and indices are reported as such in *_report */
 835 #define INDEX(i) (i)
 836 #endif
 837
 838 /* All output goes through the PRINTF macro.  */
 839 #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; }
 840
 841 /* ========================================================================== */
 842 /* === Prototypes of PRIVATE routines ======================================= */
 843 /* ========================================================================== */
 844
 845 PRIVATE Int init_rows_cols
 846 (
 847     Int n_row,
 848     Int n_col,
 849     Colamd_Row Row [],
 850     Colamd_Col Col [],
 851     Int A [],
 852     Int p [],
 853     Int stats [COLAMD_STATS]
 854 ) ;
 855
 856 PRIVATE void init_scoring
 857 (
 858     Int n_row,
 859     Int n_col,
 860     Colamd_Row Row [],
 861     Colamd_Col Col [],
 862     Int A [],
 863     Int head [],
 864     double knobs [COLAMD_KNOBS],
 865     Int *p_n_row2,
 866     Int *p_n_col2,
 867     Int *p_max_deg
 868 ) ;
 869
 870 PRIVATE Int find_ordering
 871 (
 872     Int n_row,
 873     Int n_col,
 874     Int Alen,
 875     Colamd_Row Row [],
 876     Colamd_Col Col [],
 877     Int A [],
 878     Int head [],
 879     Int n_col2,
 880     Int max_deg,
 881     Int pfree,
 882     Int aggressive
 883 ) ;
 884
 885 PRIVATE void order_children
 886 (
 887     Int n_col,
 888     Colamd_Col Col [],
 889     Int p []
 890 ) ;
 891
 892 PRIVATE void detect_super_cols
 893 (
 894
 895 #ifndef NDEBUG
 896     Int n_col,
 897     Colamd_Row Row [],
 898 #endif /* NDEBUG */
 899
 900     Colamd_Col Col [],
 901     Int A [],
 902     Int head [],
 903     Int row_start,
 904     Int row_length
 905 ) ;
 906
 907 PRIVATE Int garbage_collection
 908 (
 909     Int n_row,
 910     Int n_col,
 911     Colamd_Row Row [],
 912     Colamd_Col Col [],
 913     Int A [],
 914     Int *pfree
 915 ) ;
 916
 917 PRIVATE Int clear_mark
 918 (
 919     Int tag_mark,
 920     Int max_mark,
 921     Int n_row,
 922     Colamd_Row Row []
 923 ) ;
 924
 925 PRIVATE void print_report
 926 (
 927     char *method,
 928     Int stats [COLAMD_STATS]
 929 ) ;
 930
 931 /* ========================================================================== */
 932 /* === Debugging prototypes and definitions ================================= */
 933 /* ========================================================================== */
 934
 935 #ifndef NDEBUG
 936
 937 #include <assert.h>
 938
 939 /* colamd_debug is the *ONLY* global variable, and is only */
 940 /* present when debugging */
 941
 942 PRIVATE Int colamd_debug = 0 ;  /* debug print level */
 943
 944 #define DEBUG0(params) { PRINTF (params) ; }
 945 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; }
 946 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; }
 947 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; }
 948 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; }
 949
 950 #ifdef MATLAB_MEX_FILE
 951 #define ASSERT(expression) (mxAssert ((expression), ""))
 952 #else
 953 #define ASSERT(expression) (assert (expression))
 954 #endif /* MATLAB_MEX_FILE */
 955
 956 PRIVATE void colamd_get_debug   /* gets the debug print level from getenv */
 957 (
 958     char *method
 959 ) ;
 960
 961 PRIVATE void debug_deg_lists
 962 (
 963     Int n_row,
 964     Int n_col,
 965     Colamd_Row Row [],
 966     Colamd_Col Col [],
 967     Int head [],
 968     Int min_score,
 969     Int should,
 970     Int max_deg
 971 ) ;
 972
 973 PRIVATE void debug_mark
 974 (
 975     Int n_row,
 976     Colamd_Row Row [],
 977     Int tag_mark,
 978     Int max_mark
 979 ) ;
 980
 981 PRIVATE void debug_matrix
 982 (
 983     Int n_row,
 984     Int n_col,
 985     Colamd_Row Row [],
 986     Colamd_Col Col [],
 987     Int A []
 988 ) ;
 989
 990 PRIVATE void debug_structures
 991 (
 992     Int n_row,
 993     Int n_col,
 994     Colamd_Row Row [],
 995     Colamd_Col Col [],
 996     Int A [],
 997     Int n_col2
 998 ) ;
 999
1000 #else /* NDEBUG */
1001
1002 /* === No debugging ========================================================= */
1003
1004 #define DEBUG0(params) ;
1005 #define DEBUG1(params) ;
1006 #define DEBUG2(params) ;
1007 #define DEBUG3(params) ;
1008 #define DEBUG4(params) ;
1009
1010 #define ASSERT(expression)
1011
1012 #endif /* NDEBUG */
1013
1014 /* ========================================================================== */
1015 /* === USER-CALLABLE ROUTINES: ============================================== */
1016 /* ========================================================================== */
1017
1018 /* ========================================================================== */
1019 /* === colamd_recommended =================================================== */
1020 /* ========================================================================== */
1021
1022 /*
1023     The colamd_recommended routine returns the suggested size for Alen.  This
1024     value has been determined to provide good balance between the number of
1025     garbage collections and the memory requirements for colamd.  If any
1026     argument is negative, or if integer overflow occurs, a 0 is returned as an
1027     error condition.  2*nnz space is required for the row and column
1028     indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is
1029     required for the Col and Row arrays, respectively, which are internal to
1030     colamd (roughly 6*n_col + 4*n_row).  An additional n_col space is the
1031     minimal amount of "elbow room", and nnz/5 more space is recommended for
1032     run time efficiency.
1033
1034     Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10.
1035
1036     This function is not needed when using symamd.
1037 */
1038
1039 /* add two values of type size_t, and check for integer overflow */
1040 static size_t t_add (size_t a, size_t b, int *ok)
1041 {
1042     (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
1043     return ((*ok) ? (a + b) : 0) ;
1044 }
1045
1046 /* compute a*k where k is a small integer, and check for integer overflow */
1047 static size_t t_mult (size_t a, size_t k, int *ok)
1048 {
1049     size_t i, s = 0 ;
1050     for (i = 0 ; i < k ; i++)
1051     {
1052         s = t_add (s, a, ok) ;
1053     }
1054     return (s) ;
1055 }
1056
1057 /* size of the Col and Row structures */
1058 #define COLAMD_C(n_col,ok) \
1059     ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int)))
1060
1061 #define COLAMD_R(n_row,ok) \
1062     ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int)))
1063
1064
1065 PUBLIC size_t COLAMD_recommended        /* returns recommended value of Alen. */
1066 (
1067     /* === Parameters ======================================================= */
1068
1069     Int nnz,                    /* number of nonzeros in A */
1070     Int n_row,                  /* number of rows in A */
1071     Int n_col                   /* number of columns in A */
1072 )
1073 {
1074     size_t s, c, r ;
1075     int ok = TRUE ;
1076     if (nnz < 0 || n_row < 0 || n_col < 0)
1077     {
1078         return (0) ;
1079     }
1080     s = t_mult (nnz, 2, &ok) ;      /* 2*nnz */
1081     c = COLAMD_C (n_col, &ok) ;     /* size of column structures */
1082     r = COLAMD_R (n_row, &ok) ;     /* size of row structures */
1083     s = t_add (s, c, &ok) ;
1084     s = t_add (s, r, &ok) ;
1085     s = t_add (s, n_col, &ok) ;     /* elbow room */
1086     s = t_add (s, nnz/5, &ok) ;     /* elbow room */
1087     ok = ok && (s < Int_MAX) ;
1088     return (ok ? s : 0) ;
1089 }
1090
1091
1092 /* ========================================================================== */
1093 /* === colamd_set_defaults ================================================== */
1094 /* ========================================================================== */
1095
1096 /*
1097     The colamd_set_defaults routine sets the default values of the user-
1098     controllable parameters for colamd and symamd:
1099
1100         Colamd: rows with more than max (16, knobs [0] * sqrt (n_col))
1101         entries are removed prior to ordering.  Columns with more than
1102         max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed
1103         prior to ordering, and placed last in the output column ordering.
1104
1105         Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n))
1106         entries are removed prior to ordering, and placed last in the
1107         output ordering.
1108
1109         knobs [0]       dense row control
1110
1111         knobs [1]       dense column control
1112
1113         knobs [2]       if nonzero, do aggresive absorption
1114
1115         knobs [3..19]   unused, but future versions might use this
1116
1117 */
1118
1119 PUBLIC void COLAMD_set_defaults
1120 (
1121     /* === Parameters ======================================================= */
1122
1123     double knobs [COLAMD_KNOBS]         /* knob array */
1124 )
1125 {
1126     /* === Local variables ================================================== */
1127
1128     Int i ;
1129
1130     if (!knobs)
1131     {
1132         return ;                        /* no knobs to initialize */
1133     }
1134     for (i = 0 ; i < COLAMD_KNOBS ; i++)
1135     {
1136         knobs [i] = 0 ;
1137     }
1138     knobs [COLAMD_DENSE_ROW] = 10 ;
1139     knobs [COLAMD_DENSE_COL] = 10 ;
1140     knobs [COLAMD_AGGRESSIVE] = TRUE ;  /* default: do aggressive absorption*/
1141 }
1142
1143
1144 /* ========================================================================== */
1145 /* === symamd =============================================================== */
1146 /* ========================================================================== */
1147
1148 PUBLIC Int SYMAMD_MAIN                  /* return TRUE if OK, FALSE otherwise */
1149 (
1150     /* === Parameters ======================================================= */
1151
1152     Int n,                              /* number of rows and columns of A */
1153     Int A [],                           /* row indices of A */
1154     Int p [],                           /* column pointers of A */
1155     Int perm [],                        /* output permutation, size n+1 */
1156     double knobs [COLAMD_KNOBS],        /* parameters (uses defaults if NULL) */
1157     Int stats [COLAMD_STATS],           /* output statistics and error codes */
1158     void * (*allocate) (size_t, size_t),
1159                                         /* pointer to calloc (ANSI C) or */
1160                                         /* mxCalloc (for MATLAB mexFunction) */
1161     void (*release) (void *)
1162                                         /* pointer to free (ANSI C) or */
1163                                         /* mxFree (for MATLAB mexFunction) */
1164 )
1165 {
1166     /* === Local variables ================================================== */
1167
1168     Int *count ;                /* length of each column of M, and col pointer*/
1169     Int *mark ;                 /* mark array for finding duplicate entries */
1170     Int *M ;                    /* row indices of matrix M */
1171     size_t Mlen ;               /* length of M */
1172     Int n_row ;                 /* number of rows in M */
1173     Int nnz ;                   /* number of entries in A */
1174     Int i ;                     /* row index of A */
1175     Int j ;                     /* column index of A */
1176     Int k ;                     /* row index of M */
1177     Int mnz ;                   /* number of nonzeros in M */
1178     Int pp ;                    /* index into a column of A */
1179     Int last_row ;              /* last row seen in the current column */
1180     Int length ;                /* number of nonzeros in a column */
1181
1182     double cknobs [COLAMD_KNOBS] ;              /* knobs for colamd */
1183     double default_knobs [COLAMD_KNOBS] ;       /* default knobs for colamd */
1184
1185 #ifndef NDEBUG
1186     colamd_get_debug ("symamd") ;
1187 #endif /* NDEBUG */
1188
1189     /* === Check the input arguments ======================================== */
1190
1191     if (!stats)
1192     {
1193         DEBUG0 (("symamd: stats not present\n")) ;
1194         return (FALSE) ;
1195     }
1196     for (i = 0 ; i < COLAMD_STATS ; i++)
1197     {
1198         stats [i] = 0 ;
1199     }
1200     stats [COLAMD_STATUS] = COLAMD_OK ;
1201     stats [COLAMD_INFO1] = -1 ;
1202     stats [COLAMD_INFO2] = -1 ;
1203
1204     if (!A)
1205     {
1206         stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1207         DEBUG0 (("symamd: A not present\n")) ;
1208         return (FALSE) ;
1209     }
1210
1211     if (!p)             /* p is not present */
1212     {
1213         stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1214         DEBUG0 (("symamd: p not present\n")) ;
1215         return (FALSE) ;
1216     }
1217
1218     if (n < 0)          /* n must be >= 0 */
1219     {
1220         stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1221         stats [COLAMD_INFO1] = n ;
1222         DEBUG0 (("symamd: n negative %d\n", n)) ;
1223         return (FALSE) ;
1224     }
1225
1226     nnz = p [n] ;
1227     if (nnz < 0)        /* nnz must be >= 0 */
1228     {
1229         stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1230         stats [COLAMD_INFO1] = nnz ;
1231         DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ;
1232         return (FALSE) ;
1233     }
1234
1235     if (p [0] != 0)
1236     {
1237         stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1238         stats [COLAMD_INFO1] = p [0] ;
1239         DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ;
1240         return (FALSE) ;
1241     }
1242
1243     /* === If no knobs, set default knobs =================================== */
1244
1245     if (!knobs)
1246     {
1247         COLAMD_set_defaults (default_knobs) ;
1248         knobs = default_knobs ;
1249     }
1250
1251     /* === Allocate count and mark ========================================== */
1252
1253     count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1254     if (!count)
1255     {
1256         stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1257         DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ;
1258         return (FALSE) ;
1259     }
1260
1261     mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1262     if (!mark)
1263     {
1264         stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1265         (*release) ((void *) count) ;
1266         DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ;
1267         return (FALSE) ;
1268     }
1269
1270     /* === Compute column counts of M, check if A is valid ================== */
1271
1272     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
1273
1274     for (i = 0 ; i < n ; i++)
1275     {
1276         mark [i] = -1 ;
1277     }
1278
1279     for (j = 0 ; j < n ; j++)
1280     {
1281         last_row = -1 ;
1282
1283         length = p [j+1] - p [j] ;
1284         if (length < 0)
1285         {
1286             /* column pointers must be non-decreasing */
1287             stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1288             stats [COLAMD_INFO1] = j ;
1289             stats [COLAMD_INFO2] = length ;
1290             (*release) ((void *) count) ;
1291             (*release) ((void *) mark) ;
1292             DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ;
1293             return (FALSE) ;
1294         }
1295
1296         for (pp = p [j] ; pp < p [j+1] ; pp++)
1297         {
1298             i = A [pp] ;
1299             if (i < 0 || i >= n)
1300             {
1301                 /* row index i, in column j, is out of bounds */
1302                 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1303                 stats [COLAMD_INFO1] = j ;
1304                 stats [COLAMD_INFO2] = i ;
1305                 stats [COLAMD_INFO3] = n ;
1306                 (*release) ((void *) count) ;
1307                 (*release) ((void *) mark) ;
1308                 DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ;
1309                 return (FALSE) ;
1310             }
1311
1312             if (i <= last_row || mark [i] == j)
1313             {
1314                 /* row index is unsorted or repeated (or both), thus col */
1315                 /* is jumbled.  This is a notice, not an error condition. */
1316                 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1317                 stats [COLAMD_INFO1] = j ;
1318                 stats [COLAMD_INFO2] = i ;
1319                 (stats [COLAMD_INFO3]) ++ ;
1320                 DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
1321             }
1322
1323             if (i > j && mark [i] != j)
1324             {
1325                 /* row k of M will contain column indices i and j */
1326                 count [i]++ ;
1327                 count [j]++ ;
1328             }
1329
1330             /* mark the row as having been seen in this column */
1331             mark [i] = j ;
1332
1333             last_row = i ;
1334         }
1335     }
1336
1337     /* v2.4: removed free(mark) */
1338
1339     /* === Compute column pointers of M ===================================== */
1340
1341     /* use output permutation, perm, for column pointers of M */
1342     perm [0] = 0 ;
1343     for (j = 1 ; j <= n ; j++)
1344     {
1345         perm [j] = perm [j-1] + count [j-1] ;
1346     }
1347     for (j = 0 ; j < n ; j++)
1348     {
1349         count [j] = perm [j] ;
1350     }
1351
1352     /* === Construct M ====================================================== */
1353
1354     mnz = perm [n] ;
1355     n_row = mnz / 2 ;
1356     Mlen = COLAMD_recommended (mnz, n_row, n) ;
1357     M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
1358     DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n",
1359         n_row, n, mnz, (double) Mlen)) ;
1360
1361     if (!M)
1362     {
1363         stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ;
1364         (*release) ((void *) count) ;
1365         (*release) ((void *) mark) ;
1366         DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ;
1367         return (FALSE) ;
1368     }
1369
1370     k = 0 ;
1371
1372     if (stats [COLAMD_STATUS] == COLAMD_OK)
1373     {
1374         /* Matrix is OK */
1375         for (j = 0 ; j < n ; j++)
1376         {
1377             ASSERT (p [j+1] - p [j] >= 0) ;
1378             for (pp = p [j] ; pp < p [j+1] ; pp++)
1379             {
1380                 i = A [pp] ;
1381                 ASSERT (i >= 0 && i < n) ;
1382                 if (i > j)
1383                 {
1384                     /* row k of M contains column indices i and j */
1385                     M [count [i]++] = k ;
1386                     M [count [j]++] = k ;
1387                     k++ ;
1388                 }
1389             }
1390         }
1391     }
1392     else
1393     {
1394         /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
1395         DEBUG0 (("symamd: Duplicates in A.\n")) ;
1396         for (i = 0 ; i < n ; i++)
1397         {
1398             mark [i] = -1 ;
1399         }
1400         for (j = 0 ; j < n ; j++)
1401         {
1402             ASSERT (p [j+1] - p [j] >= 0) ;
1403             for (pp = p [j] ; pp < p [j+1] ; pp++)
1404             {
1405                 i = A [pp] ;
1406                 ASSERT (i >= 0 && i < n) ;
1407                 if (i > j && mark [i] != j)
1408                 {
1409                     /* row k of M contains column indices i and j */
1410                     M [count [i]++] = k ;
1411                     M [count [j]++] = k ;
1412                     k++ ;
1413                     mark [i] = j ;
1414                 }
1415             }
1416         }
1417         /* v2.4: free(mark) moved below */
1418     }
1419
1420     /* count and mark no longer needed */
1421     (*release) ((void *) count) ;
1422     (*release) ((void *) mark) ;        /* v2.4: free (mark) moved here */
1423     ASSERT (k == n_row) ;
1424
1425     /* === Adjust the knobs for M =========================================== */
1426
1427     for (i = 0 ; i < COLAMD_KNOBS ; i++)
1428     {
1429         cknobs [i] = knobs [i] ;
1430     }
1431
1432     /* there are no dense rows in M */
1433     cknobs [COLAMD_DENSE_ROW] = -1 ;
1434     cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ;
1435
1436     /* === Order the columns of M =========================================== */
1437
1438     /* v2.4: colamd cannot fail here, so the error check is removed */
1439     (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ;
1440
1441     /* Note that the output permutation is now in perm */
1442
1443     /* === get the statistics for symamd from colamd ======================== */
1444
1445     /* a dense column in colamd means a dense row and col in symamd */
1446     stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ;
1447
1448     /* === Free M =========================================================== */
1449
1450     (*release) ((void *) M) ;
1451     DEBUG0 (("symamd: done.\n")) ;
1452     return (TRUE) ;
1453
1454 }
1455
1456 /* ========================================================================== */
1457 /* === colamd =============================================================== */
1458 /* ========================================================================== */
1459
1460 /*
1461     The colamd routine computes a column ordering Q of a sparse matrix
1462     A such that the LU factorization P(AQ) = LU remains sparse, where P is
1463     selected via partial pivoting.   The routine can also be viewed as
1464     providing a permutation Q such that the Cholesky factorization
1465     (AQ)'(AQ) = LL' remains sparse.
1466 */
1467
1468 PUBLIC Int COLAMD_MAIN          /* returns TRUE if successful, FALSE otherwise*/
1469 (
1470     /* === Parameters ======================================================= */
1471
1472     Int n_row,                  /* number of rows in A */
1473     Int n_col,                  /* number of columns in A */
1474     Int Alen,                   /* length of A */
1475     Int A [],                   /* row indices of A */
1476     Int p [],                   /* pointers to columns in A */
1477     double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */
1478     Int stats [COLAMD_STATS]    /* output statistics and error codes */
1479 )
1480 {
1481     /* === Local variables ================================================== */
1482
1483     Int i ;                     /* loop index */
1484     Int nnz ;                   /* nonzeros in A */
1485     size_t Row_size ;           /* size of Row [], in integers */
1486     size_t Col_size ;           /* size of Col [], in integers */
1487     size_t need ;               /* minimum required length of A */
1488     Colamd_Row *Row ;           /* pointer into A of Row [0..n_row] array */
1489     Colamd_Col *Col ;           /* pointer into A of Col [0..n_col] array */
1490     Int n_col2 ;                /* number of non-dense, non-empty columns */
1491     Int n_row2 ;                /* number of non-dense, non-empty rows */
1492     Int ngarbage ;              /* number of garbage collections performed */
1493     Int max_deg ;               /* maximum row degree */
1494     double default_knobs [COLAMD_KNOBS] ;       /* default knobs array */
1495     Int aggressive ;            /* do aggressive absorption */
1496     int ok ;
1497
1498 #ifndef NDEBUG
1499     colamd_get_debug ("colamd") ;
1500 #endif /* NDEBUG */
1501
1502     /* === Check the input arguments ======================================== */
1503
1504     if (!stats)
1505     {
1506         DEBUG0 (("colamd: stats not present\n")) ;
1507         return (FALSE) ;
1508     }
1509     for (i = 0 ; i < COLAMD_STATS ; i++)
1510     {
1511         stats [i] = 0 ;
1512     }
1513     stats [COLAMD_STATUS] = COLAMD_OK ;
1514     stats [COLAMD_INFO1] = -1 ;
1515     stats [COLAMD_INFO2] = -1 ;
1516
1517     if (!A)             /* A is not present */
1518     {
1519         stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
1520         DEBUG0 (("colamd: A not present\n")) ;
1521         return (FALSE) ;
1522     }
1523
1524     if (!p)             /* p is not present */
1525     {
1526         stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
1527         DEBUG0 (("colamd: p not present\n")) ;
1528         return (FALSE) ;
1529     }
1530
1531     if (n_row < 0)      /* n_row must be >= 0 */
1532     {
1533         stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
1534         stats [COLAMD_INFO1] = n_row ;
1535         DEBUG0 (("colamd: nrow negative %d\n", n_row)) ;
1536         return (FALSE) ;
1537     }
1538
1539     if (n_col < 0)      /* n_col must be >= 0 */
1540     {
1541         stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
1542         stats [COLAMD_INFO1] = n_col ;
1543         DEBUG0 (("colamd: ncol negative %d\n", n_col)) ;
1544         return (FALSE) ;
1545     }
1546
1547     nnz = p [n_col] ;
1548     if (nnz < 0)        /* nnz must be >= 0 */
1549     {
1550         stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
1551         stats [COLAMD_INFO1] = nnz ;
1552         DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ;
1553         return (FALSE) ;
1554     }
1555
1556     if (p [0] != 0)
1557     {
1558         stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
1559         stats [COLAMD_INFO1] = p [0] ;
1560         DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ;
1561         return (FALSE) ;
1562     }
1563
1564     /* === If no knobs, set default knobs =================================== */
1565
1566     if (!knobs)
1567     {
1568         COLAMD_set_defaults (default_knobs) ;
1569         knobs = default_knobs ;
1570     }
1571
1572     aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ;
1573
1574     /* === Allocate the Row and Col arrays from array A ===================== */
1575
1576     ok = TRUE ;
1577     Col_size = COLAMD_C (n_col, &ok) ;      /* size of Col array of structs */
1578     Row_size = COLAMD_R (n_row, &ok) ;      /* size of Row array of structs */
1579
1580     /* need = 2*nnz + n_col + Col_size + Row_size ; */
1581     need = t_mult (nnz, 2, &ok) ;
1582     need = t_add (need, n_col, &ok) ;
1583     need = t_add (need, Col_size, &ok) ;
1584     need = t_add (need, Row_size, &ok) ;
1585
1586     if (!ok || need > (size_t) Alen || need > Int_MAX)
1587     {
1588         /* not enough space in array A to perform the ordering */
1589         stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
1590         stats [COLAMD_INFO1] = need ;
1591         stats [COLAMD_INFO2] = Alen ;
1592         DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
1593         return (FALSE) ;
1594     }
1595
1596     Alen -= Col_size + Row_size ;
1597     Col = (Colamd_Col *) &A [Alen] ;
1598     Row = (Colamd_Row *) &A [Alen + Col_size] ;
1599
1600     /* === Construct the row and column data structures ===================== */
1601
1602     if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
1603     {
1604         /* input matrix is invalid */
1605         DEBUG0 (("colamd: Matrix invalid\n")) ;
1606         return (FALSE) ;
1607     }
1608
1609     /* === Initialize scores, kill dense rows/columns ======================= */
1610
1611     init_scoring (n_row, n_col, Row, Col, A, p, knobs,
1612         &n_row2, &n_col2, &max_deg) ;
1613
1614     /* === Order the supercolumns =========================================== */
1615
1616     ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
1617         n_col2, max_deg, 2*nnz, aggressive) ;
1618
1619     /* === Order the non-principal columns ================================== */
1620
1621     order_children (n_col, Col, p) ;
1622
1623     /* === Return statistics in stats ======================================= */
1624
1625     stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
1626     stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
1627     stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
1628     DEBUG0 (("colamd: done.\n")) ;
1629     return (TRUE) ;
1630 }
1631
1632
1633 /* ========================================================================== */
1634 /* === colamd_report ======================================================== */
1635 /* ========================================================================== */
1636
1637 PUBLIC void COLAMD_report
1638 (
1639     Int stats [COLAMD_STATS]
1640 )
1641 {
1642     print_report ("colamd", stats) ;
1643 }
1644
1645
1646 /* ========================================================================== */
1647 /* === symamd_report ======================================================== */
1648 /* ========================================================================== */
1649
1650 PUBLIC void SYMAMD_report
1651 (
1652     Int stats [COLAMD_STATS]
1653 )
1654 {
1655     print_report ("symamd", stats) ;
1656 }
1657
1658
1659
1660 /* ========================================================================== */
1661 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
1662 /* ========================================================================== */
1663
1664 /* There are no user-callable routines beyond this point in the file */
1665
1666
1667 /* ========================================================================== */
1668 /* === init_rows_cols ======================================================= */
1669 /* ========================================================================== */
1670
1671 /*
1672     Takes the column form of the matrix in A and creates the row form of the
1673     matrix.  Also, row and column attributes are stored in the Col and Row
1674     structs.  If the columns are un-sorted or contain duplicate row indices,
1675     this routine will also sort and remove duplicate row indices from the
1676     column form of the matrix.  Returns FALSE if the matrix is invalid,
1677     TRUE otherwise.  Not user-callable.
1678 */
1679
1680 PRIVATE Int init_rows_cols      /* returns TRUE if OK, or FALSE otherwise */
1681 (
1682     /* === Parameters ======================================================= */
1683
1684     Int n_row,                  /* number of rows of A */
1685     Int n_col,                  /* number of columns of A */
1686     Colamd_Row Row [],          /* of size n_row+1 */
1687     Colamd_Col Col [],          /* of size n_col+1 */
1688     Int A [],                   /* row indices of A, of size Alen */
1689     Int p [],                   /* pointers to columns in A, of size n_col+1 */
1690     Int stats [COLAMD_STATS]    /* colamd statistics */
1691 )
1692 {
1693     /* === Local variables ================================================== */
1694
1695     Int col ;                   /* a column index */
1696     Int row ;                   /* a row index */
1697     Int *cp ;                   /* a column pointer */
1698     Int *cp_end ;               /* a pointer to the end of a column */
1699     Int *rp ;                   /* a row pointer */
1700     Int *rp_end ;               /* a pointer to the end of a row */
1701     Int last_row ;              /* previous row */
1702
1703     /* === Initialize columns, and check column pointers ==================== */
1704
1705     for (col = 0 ; col < n_col ; col++)
1706     {
1707         Col [col].start = p [col] ;
1708         Col [col].length = p [col+1] - p [col] ;
1709
1710         if (Col [col].length < 0)
1711         {
1712             /* column pointers must be non-decreasing */
1713             stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
1714             stats [COLAMD_INFO1] = col ;
1715             stats [COLAMD_INFO2] = Col [col].length ;
1716             DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ;
1717             return (FALSE) ;
1718         }
1719
1720         Col [col].shared1.thickness = 1 ;
1721         Col [col].shared2.score = 0 ;
1722         Col [col].shared3.prev = EMPTY ;
1723         Col [col].shared4.degree_next = EMPTY ;
1724     }
1725
1726     /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
1727
1728     /* === Scan columns, compute row degrees, and check row indices ========= */
1729
1730     stats [COLAMD_INFO3] = 0 ;  /* number of duplicate or unsorted row indices*/
1731
1732     for (row = 0 ; row < n_row ; row++)
1733     {
1734         Row [row].length = 0 ;
1735         Row [row].shared2.mark = -1 ;
1736     }
1737
1738     for (col = 0 ; col < n_col ; col++)
1739     {
1740         last_row = -1 ;
1741
1742         cp = &A [p [col]] ;
1743         cp_end = &A [p [col+1]] ;
1744
1745         while (cp < cp_end)
1746         {
1747             row = *cp++ ;
1748
1749             /* make sure row indices within range */
1750             if (row < 0 || row >= n_row)
1751             {
1752                 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
1753                 stats [COLAMD_INFO1] = col ;
1754                 stats [COLAMD_INFO2] = row ;
1755                 stats [COLAMD_INFO3] = n_row ;
1756                 DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ;
1757                 return (FALSE) ;
1758             }
1759
1760             if (row <= last_row || Row [row].shared2.mark == col)
1761             {
1762                 /* row index are unsorted or repeated (or both), thus col */
1763                 /* is jumbled.  This is a notice, not an error condition. */
1764                 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
1765                 stats [COLAMD_INFO1] = col ;
1766                 stats [COLAMD_INFO2] = row ;
1767                 (stats [COLAMD_INFO3]) ++ ;
1768                 DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col));
1769             }
1770
1771             if (Row [row].shared2.mark != col)
1772             {
1773                 Row [row].length++ ;
1774             }
1775             else
1776             {
1777                 /* this is a repeated entry in the column, */
1778                 /* it will be removed */
1779                 Col [col].length-- ;
1780             }
1781
1782             /* mark the row as having been seen in this column */
1783             Row [row].shared2.mark = col ;
1784
1785             last_row = row ;
1786         }
1787     }
1788
1789     /* === Compute row pointers ============================================= */
1790
1791     /* row form of the matrix starts directly after the column */
1792     /* form of matrix in A */
1793     Row [0].start = p [n_col] ;
1794     Row [0].shared1.p = Row [0].start ;
1795     Row [0].shared2.mark = -1 ;
1796     for (row = 1 ; row < n_row ; row++)
1797     {
1798         Row [row].start = Row [row-1].start + Row [row-1].length ;
1799         Row [row].shared1.p = Row [row].start ;
1800         Row [row].shared2.mark = -1 ;
1801     }
1802
1803     /* === Create row form ================================================== */
1804
1805     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1806     {
1807         /* if cols jumbled, watch for repeated row indices */
1808         for (col = 0 ; col < n_col ; col++)
1809         {
1810             cp = &A [p [col]] ;
1811             cp_end = &A [p [col+1]] ;
1812             while (cp < cp_end)
1813             {
1814                 row = *cp++ ;
1815                 if (Row [row].shared2.mark != col)
1816                 {
1817                     A [(Row [row].shared1.p)++] = col ;
1818                     Row [row].shared2.mark = col ;
1819                 }
1820             }
1821         }
1822     }
1823     else
1824     {
1825         /* if cols not jumbled, we don't need the mark (this is faster) */
1826         for (col = 0 ; col < n_col ; col++)
1827         {
1828             cp = &A [p [col]] ;
1829             cp_end = &A [p [col+1]] ;
1830             while (cp < cp_end)
1831             {
1832                 A [(Row [*cp++].shared1.p)++] = col ;
1833             }
1834         }
1835     }
1836
1837     /* === Clear the row marks and set row degrees ========================== */
1838
1839     for (row = 0 ; row < n_row ; row++)
1840     {
1841         Row [row].shared2.mark = 0 ;
1842         Row [row].shared1.degree = Row [row].length ;
1843     }
1844
1845     /* === See if we need to re-create columns ============================== */
1846
1847     if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
1848     {
1849         DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ;
1850
1851 #ifndef NDEBUG
1852         /* make sure column lengths are correct */
1853         for (col = 0 ; col < n_col ; col++)
1854         {
1855             p [col] = Col [col].length ;
1856         }
1857         for (row = 0 ; row < n_row ; row++)
1858         {
1859             rp = &A [Row [row].start] ;
1860             rp_end = rp + Row [row].length ;
1861             while (rp < rp_end)
1862             {
1863                 p [*rp++]-- ;
1864             }
1865         }
1866         for (col = 0 ; col < n_col ; col++)
1867         {
1868             ASSERT (p [col] == 0) ;
1869         }
1870         /* now p is all zero (different than when debugging is turned off) */
1871 #endif /* NDEBUG */
1872
1873         /* === Compute col pointers ========================================= */
1874
1875         /* col form of the matrix starts at A [0]. */
1876         /* Note, we may have a gap between the col form and the row */
1877         /* form if there were duplicate entries, if so, it will be */
1878         /* removed upon the first garbage collection */
1879         Col [0].start = 0 ;
1880         p [0] = Col [0].start ;
1881         for (col = 1 ; col < n_col ; col++)
1882         {
1883             /* note that the lengths here are for pruned columns, i.e. */
1884             /* no duplicate row indices will exist for these columns */
1885             Col [col].start = Col [col-1].start + Col [col-1].length ;
1886             p [col] = Col [col].start ;
1887         }
1888
1889         /* === Re-create col form =========================================== */
1890
1891         for (row = 0 ; row < n_row ; row++)
1892         {
1893             rp = &A [Row [row].start] ;
1894             rp_end = rp + Row [row].length ;
1895             while (rp < rp_end)
1896             {
1897                 A [(p [*rp++])++] = row ;
1898             }
1899         }
1900     }
1901
1902     /* === Done.  Matrix is not (or no longer) jumbled ====================== */
1903
1904     return (TRUE) ;
1905 }
1906
1907
1908 /* ========================================================================== */
1909 /* === init_scoring ========================================================= */
1910 /* ========================================================================== */
1911
1912 /*
1913     Kills dense or empty columns and rows, calculates an initial score for
1914     each column, and places all columns in the degree lists.  Not user-callable.
1915 */
1916
1917 PRIVATE void init_scoring
1918 (
1919     /* === Parameters ======================================================= */
1920
1921     Int n_row,                  /* number of rows of A */
1922     Int n_col,                  /* number of columns of A */
1923     Colamd_Row Row [],          /* of size n_row+1 */
1924     Colamd_Col Col [],          /* of size n_col+1 */
1925     Int A [],                   /* column form and row form of A */
1926     Int head [],                /* of size n_col+1 */
1927     double knobs [COLAMD_KNOBS],/* parameters */
1928     Int *p_n_row2,              /* number of non-dense, non-empty rows */
1929     Int *p_n_col2,              /* number of non-dense, non-empty columns */
1930     Int *p_max_deg              /* maximum row degree */
1931 )
1932 {
1933     /* === Local variables ================================================== */
1934
1935     Int c ;                     /* a column index */
1936     Int r, row ;                /* a row index */
1937     Int *cp ;                   /* a column pointer */
1938     Int deg ;                   /* degree of a row or column */
1939     Int *cp_end ;               /* a pointer to the end of a column */
1940     Int *new_cp ;               /* new column pointer */
1941     Int col_length ;            /* length of pruned column */
1942     Int score ;                 /* current column score */
1943     Int n_col2 ;                /* number of non-dense, non-empty columns */
1944     Int n_row2 ;                /* number of non-dense, non-empty rows */
1945     Int dense_row_count ;       /* remove rows with more entries than this */
1946     Int dense_col_count ;       /* remove cols with more entries than this */
1947     Int min_score ;             /* smallest column score */
1948     Int max_deg ;               /* maximum row degree */
1949     Int next_col ;              /* Used to add to degree list.*/
1950
1951 #ifndef NDEBUG
1952     Int debug_count ;           /* debug only. */
1953 #endif /* NDEBUG */
1954
1955     /* === Extract knobs ==================================================== */
1956
1957     /* Note: if knobs contains a NaN, this is undefined: */
1958     if (knobs [COLAMD_DENSE_ROW] < 0)
1959     {
1960         /* only remove completely dense rows */
1961         dense_row_count = n_col-1 ;
1962     }
1963     else
1964     {
1965         dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ;
1966     }
1967     if (knobs [COLAMD_DENSE_COL] < 0)
1968     {
1969         /* only remove completely dense columns */
1970         dense_col_count = n_row-1 ;
1971     }
1972     else
1973     {
1974         dense_col_count =
1975             DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ;
1976     }
1977
1978     DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
1979     max_deg = 0 ;
1980     n_col2 = n_col ;
1981     n_row2 = n_row ;
1982
1983     /* === Kill empty columns =============================================== */
1984
1985     /* Put the empty columns at the end in their natural order, so that LU */
1986     /* factorization can proceed as far as possible. */
1987     for (c = n_col-1 ; c >= 0 ; c--)
1988     {
1989         deg = Col [c].length ;
1990         if (deg == 0)
1991         {
1992             /* this is a empty column, kill and order it last */
1993             Col [c].shared2.order = --n_col2 ;
1994             KILL_PRINCIPAL_COL (c) ;
1995         }
1996     }
1997     DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ;
1998
1999     /* === Kill dense columns =============================================== */
2000
2001     /* Put the dense columns at the end, in their natural order */
2002     for (c = n_col-1 ; c >= 0 ; c--)
2003     {
2004         /* skip any dead columns */
2005         if (COL_IS_DEAD (c))
2006         {
2007             continue ;
2008         }
2009         deg = Col [c].length ;
2010         if (deg > dense_col_count)
2011         {
2012             /* this is a dense column, kill and order it last */
2013             Col [c].shared2.order = --n_col2 ;
2014             /* decrement the row degrees */
2015             cp = &A [Col [c].start] ;
2016             cp_end = cp + Col [c].length ;
2017             while (cp < cp_end)
2018             {
2019                 Row [*cp++].shared1.degree-- ;
2020             }
2021             KILL_PRINCIPAL_COL (c) ;
2022         }
2023     }
2024     DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
2025
2026     /* === Kill dense and empty rows ======================================== */
2027
2028     for (r = 0 ; r < n_row ; r++)
2029     {
2030         deg = Row [r].shared1.degree ;
2031         ASSERT (deg >= 0 && deg <= n_col) ;
2032         if (deg > dense_row_count || deg == 0)
2033         {
2034             /* kill a dense or empty row */
2035             KILL_ROW (r) ;
2036             --n_row2 ;
2037         }
2038         else
2039         {
2040             /* keep track of max degree of remaining rows */
2041             max_deg = MAX (max_deg, deg) ;
2042         }
2043     }
2044     DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
2045
2046     /* === Compute initial column scores ==================================== */
2047
2048     /* At this point the row degrees are accurate.  They reflect the number */
2049     /* of "live" (non-dense) columns in each row.  No empty rows exist. */
2050     /* Some "live" columns may contain only dead rows, however.  These are */
2051     /* pruned in the code below. */
2052
2053     /* now find the initial matlab score for each column */
2054     for (c = n_col-1 ; c >= 0 ; c--)
2055     {
2056         /* skip dead column */
2057         if (COL_IS_DEAD (c))
2058         {
2059             continue ;
2060         }
2061         score = 0 ;
2062         cp = &A [Col [c].start] ;
2063         new_cp = cp ;
2064         cp_end = cp + Col [c].length ;
2065         while (cp < cp_end)
2066         {
2067             /* get a row */
2068             row = *cp++ ;
2069             /* skip if dead */
2070             if (ROW_IS_DEAD (row))
2071             {
2072                 continue ;
2073             }
2074             /* compact the column */
2075             *new_cp++ = row ;
2076             /* add row's external degree */
2077             score += Row [row].shared1.degree - 1 ;
2078             /* guard against integer overflow */
2079             score = MIN (score, n_col) ;
2080         }
2081         /* determine pruned column length */
2082         col_length = (Int) (new_cp - &A [Col [c].start]) ;
2083         if (col_length == 0)
2084         {
2085             /* a newly-made null column (all rows in this col are "dense" */
2086             /* and have already been killed) */
2087             DEBUG2 (("Newly null killed: %d\n", c)) ;
2088             Col [c].shared2.order = --n_col2 ;
2089             KILL_PRINCIPAL_COL (c) ;
2090         }
2091         else
2092         {
2093             /* set column length and set score */
2094             ASSERT (score >= 0) ;
2095             ASSERT (score <= n_col) ;
2096             Col [c].length = col_length ;
2097             Col [c].shared2.score = score ;
2098         }
2099     }
2100     DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n",
2101         n_col-n_col2)) ;
2102
2103     /* At this point, all empty rows and columns are dead.  All live columns */
2104     /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
2105     /* yet).  Rows may contain dead columns, but all live rows contain at */
2106     /* least one live column. */
2107
2108 #ifndef NDEBUG
2109     debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
2110 #endif /* NDEBUG */
2111
2112     /* === Initialize degree lists ========================================== */
2113
2114 #ifndef NDEBUG
2115     debug_count = 0 ;
2116 #endif /* NDEBUG */
2117
2118     /* clear the hash buckets */
2119     for (c = 0 ; c <= n_col ; c++)
2120     {
2121         head [c] = EMPTY ;
2122     }
2123     min_score = n_col ;
2124     /* place in reverse order, so low column indices are at the front */
2125     /* of the lists.  This is to encourage natural tie-breaking */
2126     for (c = n_col-1 ; c >= 0 ; c--)
2127     {
2128         /* only add principal columns to degree lists */
2129         if (COL_IS_ALIVE (c))
2130         {
2131             DEBUG4 (("place %d score %d minscore %d ncol %d\n",
2132                 c, Col [c].shared2.score, min_score, n_col)) ;
2133
2134             /* === Add columns score to DList =============================== */
2135
2136             score = Col [c].shared2.score ;
2137
2138             ASSERT (min_score >= 0) ;
2139             ASSERT (min_score <= n_col) ;
2140             ASSERT (score >= 0) ;
2141             ASSERT (score <= n_col) ;
2142             ASSERT (head [score] >= EMPTY) ;
2143
2144             /* now add this column to dList at proper score location */
2145             next_col = head [score] ;
2146             Col [c].shared3.prev = EMPTY ;
2147             Col [c].shared4.degree_next = next_col ;
2148
2149             /* if there already was a column with the same score, set its */
2150             /* previous pointer to this new column */
2151             if (next_col != EMPTY)
2152             {
2153                 Col [next_col].shared3.prev = c ;
2154             }
2155             head [score] = c ;
2156
2157             /* see if this score is less than current min */
2158             min_score = MIN (min_score, score) ;
2159
2160 #ifndef NDEBUG
2161             debug_count++ ;
2162 #endif /* NDEBUG */
2163
2164         }
2165     }
2166
2167 #ifndef NDEBUG
2168     DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n",
2169         debug_count, n_col, n_col-debug_count)) ;
2170     ASSERT (debug_count == n_col2) ;
2171     debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
2172 #endif /* NDEBUG */
2173
2174     /* === Return number of remaining columns, and max row degree =========== */
2175
2176     *p_n_col2 = n_col2 ;
2177     *p_n_row2 = n_row2 ;
2178     *p_max_deg = max_deg ;
2179 }
2180
2181
2182 /* ========================================================================== */
2183 /* === find_ordering ======================================================== */
2184 /* ========================================================================== */
2185
2186 /*
2187     Order the principal columns of the supercolumn form of the matrix
2188     (no supercolumns on input).  Uses a minimum approximate column minimum
2189     degree ordering method.  Not user-callable.
2190 */
2191
2192 PRIVATE Int find_ordering       /* return the number of garbage collections */
2193 (
2194     /* === Parameters ======================================================= */
2195
2196     Int n_row,                  /* number of rows of A */
2197     Int n_col,                  /* number of columns of A */
2198     Int Alen,                   /* size of A, 2*nnz + n_col or larger */
2199     Colamd_Row Row [],          /* of size n_row+1 */
2200     Colamd_Col Col [],          /* of size n_col+1 */
2201     Int A [],                   /* column form and row form of A */
2202     Int head [],                /* of size n_col+1 */
2203     Int n_col2,                 /* Remaining columns to order */
2204     Int max_deg,                /* Maximum row degree */
2205     Int pfree,                  /* index of first free slot (2*nnz on entry) */
2206     Int aggressive
2207 )
2208 {
2209     /* === Local variables ================================================== */
2210
2211     Int k ;                     /* current pivot ordering step */
2212     Int pivot_col ;             /* current pivot column */
2213     Int *cp ;                   /* a column pointer */
2214     Int *rp ;                   /* a row pointer */
2215     Int pivot_row ;             /* current pivot row */
2216     Int *new_cp ;               /* modified column pointer */
2217     Int *new_rp ;               /* modified row pointer */
2218     Int pivot_row_start ;       /* pointer to start of pivot row */
2219     Int pivot_row_degree ;      /* number of columns in pivot row */
2220     Int pivot_row_length ;      /* number of supercolumns in pivot row */
2221     Int pivot_col_score ;       /* score of pivot column */
2222     Int needed_memory ;         /* free space needed for pivot row */
2223     Int *cp_end ;               /* pointer to the end of a column */
2224     Int *rp_end ;               /* pointer to the end of a row */
2225     Int row ;                   /* a row index */
2226     Int col ;                   /* a column index */
2227     Int max_score ;             /* maximum possible score */
2228     Int cur_score ;             /* score of current column */
2229     unsigned Int hash ;         /* hash value for supernode detection */
2230     Int head_column ;           /* head of hash bucket */
2231     Int first_col ;             /* first column in hash bucket */
2232     Int tag_mark ;              /* marker value for mark array */
2233     Int row_mark ;              /* Row [row].shared2.mark */
2234     Int set_difference ;        /* set difference size of row with pivot row */
2235     Int min_score ;             /* smallest column score */
2236     Int col_thickness ;         /* "thickness" (no. of columns in a supercol) */
2237     Int max_mark ;              /* maximum value of tag_mark */
2238     Int pivot_col_thickness ;   /* number of columns represented by pivot col */
2239     Int prev_col ;              /* Used by Dlist operations. */
2240     Int next_col ;              /* Used by Dlist operations. */
2241     Int ngarbage ;              /* number of garbage collections performed */
2242
2243 #ifndef NDEBUG
2244     Int debug_d ;               /* debug loop counter */
2245     Int debug_step = 0 ;        /* debug loop counter */
2246 #endif /* NDEBUG */
2247
2248     /* === Initialization and clear mark ==================================== */
2249
2250     max_mark = INT_MAX - n_col ;        /* INT_MAX defined in <limits.h> */
2251     tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2252     min_score = 0 ;
2253     ngarbage = 0 ;
2254     DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ;
2255
2256     /* === Order the columns ================================================ */
2257
2258     for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */)
2259     {
2260
2261 #ifndef NDEBUG
2262         if (debug_step % 100 == 0)
2263         {
2264             DEBUG2 (("\n...       Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2265         }
2266         else
2267         {
2268             DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2269         }
2270         debug_step++ ;
2271         debug_deg_lists (n_row, n_col, Row, Col, head,
2272                 min_score, n_col2-k, max_deg) ;
2273         debug_matrix (n_row, n_col, Row, Col, A) ;
2274 #endif /* NDEBUG */
2275
2276         /* === Select pivot column, and order it ============================ */
2277
2278         /* make sure degree list isn't empty */
2279         ASSERT (min_score >= 0) ;
2280         ASSERT (min_score <= n_col) ;
2281         ASSERT (head [min_score] >= EMPTY) ;
2282
2283 #ifndef NDEBUG
2284         for (debug_d = 0 ; debug_d < min_score ; debug_d++)
2285         {
2286             ASSERT (head [debug_d] == EMPTY) ;
2287         }
2288 #endif /* NDEBUG */
2289
2290         /* get pivot column from head of minimum degree list */
2291         while (head [min_score] == EMPTY && min_score < n_col)
2292         {
2293             min_score++ ;
2294         }
2295         pivot_col = head [min_score] ;
2296         ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
2297         next_col = Col [pivot_col].shared4.degree_next ;
2298         head [min_score] = next_col ;
2299         if (next_col != EMPTY)
2300         {
2301             Col [next_col].shared3.prev = EMPTY ;
2302         }
2303
2304         ASSERT (COL_IS_ALIVE (pivot_col)) ;
2305
2306         /* remember score for defrag check */
2307         pivot_col_score = Col [pivot_col].shared2.score ;
2308
2309         /* the pivot column is the kth column in the pivot order */
2310         Col [pivot_col].shared2.order = k ;
2311
2312         /* increment order count by column thickness */
2313         pivot_col_thickness = Col [pivot_col].shared1.thickness ;
2314         k += pivot_col_thickness ;
2315         ASSERT (pivot_col_thickness > 0) ;
2316         DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ;
2317
2318         /* === Garbage_collection, if necessary ============================= */
2319
2320         needed_memory = MIN (pivot_col_score, n_col - k) ;
2321         if (pfree + needed_memory >= Alen)
2322         {
2323             pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
2324             ngarbage++ ;
2325             /* after garbage collection we will have enough */
2326             ASSERT (pfree + needed_memory < Alen) ;
2327             /* garbage collection has wiped out the Row[].shared2.mark array */
2328             tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2329
2330 #ifndef NDEBUG
2331             debug_matrix (n_row, n_col, Row, Col, A) ;
2332 #endif /* NDEBUG */
2333         }
2334
2335         /* === Compute pivot row pattern ==================================== */
2336
2337         /* get starting location for this new merged row */
2338         pivot_row_start = pfree ;
2339
2340         /* initialize new row counts to zero */
2341         pivot_row_degree = 0 ;
2342
2343         /* tag pivot column as having been visited so it isn't included */
2344         /* in merged pivot row */
2345         Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
2346
2347         /* pivot row is the union of all rows in the pivot column pattern */
2348         cp = &A [Col [pivot_col].start] ;
2349         cp_end = cp + Col [pivot_col].length ;
2350         while (cp < cp_end)
2351         {
2352             /* get a row */
2353             row = *cp++ ;
2354             DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
2355             /* skip if row is dead */
2356             if (ROW_IS_ALIVE (row))
2357             {
2358                 rp = &A [Row [row].start] ;
2359                 rp_end = rp + Row [row].length ;
2360                 while (rp < rp_end)
2361                 {
2362                     /* get a column */
2363                     col = *rp++ ;
2364                     /* add the column, if alive and untagged */
2365                     col_thickness = Col [col].shared1.thickness ;
2366                     if (col_thickness > 0 && COL_IS_ALIVE (col))
2367                     {
2368                         /* tag column in pivot row */
2369                         Col [col].shared1.thickness = -col_thickness ;
2370                         ASSERT (pfree < Alen) ;
2371                         /* place column in pivot row */
2372                         A [pfree++] = col ;
2373                         pivot_row_degree += col_thickness ;
2374                     }
2375                 }
2376             }
2377         }
2378
2379         /* clear tag on pivot column */
2380         Col [pivot_col].shared1.thickness = pivot_col_thickness ;
2381         max_deg = MAX (max_deg, pivot_row_degree) ;
2382
2383 #ifndef NDEBUG
2384         DEBUG3 (("check2\n")) ;
2385         debug_mark (n_row, Row, tag_mark, max_mark) ;
2386 #endif /* NDEBUG */
2387
2388         /* === Kill all rows used to construct pivot row ==================== */
2389
2390         /* also kill pivot row, temporarily */
2391         cp = &A [Col [pivot_col].start] ;
2392         cp_end = cp + Col [pivot_col].length ;
2393         while (cp < cp_end)
2394         {
2395             /* may be killing an already dead row */
2396             row = *cp++ ;
2397             DEBUG3 (("Kill row in pivot col: %d\n", row)) ;
2398             KILL_ROW (row) ;
2399         }
2400
2401         /* === Select a row index to use as the new pivot row =============== */
2402
2403         pivot_row_length = pfree - pivot_row_start ;
2404         if (pivot_row_length > 0)
2405         {
2406             /* pick the "pivot" row arbitrarily (first row in col) */
2407             pivot_row = A [Col [pivot_col].start] ;
2408             DEBUG3 (("Pivotal row is %d\n", pivot_row)) ;
2409         }
2410         else
2411         {
2412             /* there is no pivot row, since it is of zero length */
2413             pivot_row = EMPTY ;
2414             ASSERT (pivot_row_length == 0) ;
2415         }
2416         ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
2417
2418         /* === Approximate degree computation =============================== */
2419
2420         /* Here begins the computation of the approximate degree.  The column */
2421         /* score is the sum of the pivot row "length", plus the size of the */
2422         /* set differences of each row in the column minus the pattern of the */
2423         /* pivot row itself.  The column ("thickness") itself is also */
2424         /* excluded from the column score (we thus use an approximate */
2425         /* external degree). */
2426
2427         /* The time taken by the following code (compute set differences, and */
2428         /* add them up) is proportional to the size of the data structure */
2429         /* being scanned - that is, the sum of the sizes of each column in */
2430         /* the pivot row.  Thus, the amortized time to compute a column score */
2431         /* is proportional to the size of that column (where size, in this */
2432         /* context, is the column "length", or the number of row indices */
2433         /* in that column).  The number of row indices in a column is */
2434         /* monotonically non-decreasing, from the length of the original */
2435         /* column on input to colamd. */
2436
2437         /* === Compute set differences ====================================== */
2438
2439         DEBUG3 (("** Computing set differences phase. **\n")) ;
2440
2441         /* pivot row is currently dead - it will be revived later. */
2442
2443         DEBUG3 (("Pivot row: ")) ;
2444         /* for each column in pivot row */
2445         rp = &A [pivot_row_start] ;
2446         rp_end = rp + pivot_row_length ;
2447         while (rp < rp_end)
2448         {
2449             col = *rp++ ;
2450             ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2451             DEBUG3 (("Col: %d\n", col)) ;
2452
2453             /* clear tags used to construct pivot row pattern */
2454             col_thickness = -Col [col].shared1.thickness ;
2455             ASSERT (col_thickness > 0) ;
2456             Col [col].shared1.thickness = col_thickness ;
2457
2458             /* === Remove column from degree list =========================== */
2459
2460             cur_score = Col [col].shared2.score ;
2461             prev_col = Col [col].shared3.prev ;
2462             next_col = Col [col].shared4.degree_next ;
2463             ASSERT (cur_score >= 0) ;
2464             ASSERT (cur_score <= n_col) ;
2465             ASSERT (cur_score >= EMPTY) ;
2466             if (prev_col == EMPTY)
2467             {
2468                 head [cur_score] = next_col ;
2469             }
2470             else
2471             {
2472                 Col [prev_col].shared4.degree_next = next_col ;
2473             }
2474             if (next_col != EMPTY)
2475             {
2476                 Col [next_col].shared3.prev = prev_col ;
2477             }
2478
2479             /* === Scan the column ========================================== */
2480
2481             cp = &A [Col [col].start] ;
2482             cp_end = cp + Col [col].length ;
2483             while (cp < cp_end)
2484             {
2485                 /* get a row */
2486                 row = *cp++ ;
2487                 row_mark = Row [row].shared2.mark ;
2488                 /* skip if dead */
2489                 if (ROW_IS_MARKED_DEAD (row_mark))
2490                 {
2491                     continue ;
2492                 }
2493                 ASSERT (row != pivot_row) ;
2494                 set_difference = row_mark - tag_mark ;
2495                 /* check if the row has been seen yet */
2496                 if (set_difference < 0)
2497                 {
2498                     ASSERT (Row [row].shared1.degree <= max_deg) ;
2499                     set_difference = Row [row].shared1.degree ;
2500                 }
2501                 /* subtract column thickness from this row's set difference */
2502                 set_difference -= col_thickness ;
2503                 ASSERT (set_difference >= 0) ;
2504                 /* absorb this row if the set difference becomes zero */
2505                 if (set_difference == 0 && aggressive)
2506                 {
2507                     DEBUG3 (("aggressive absorption. Row: %d\n", row)) ;
2508                     KILL_ROW (row) ;
2509                 }
2510                 else
2511                 {
2512                     /* save the new mark */
2513                     Row [row].shared2.mark = set_difference + tag_mark ;
2514                 }
2515             }
2516         }
2517
2518 #ifndef NDEBUG
2519         debug_deg_lists (n_row, n_col, Row, Col, head,
2520                 min_score, n_col2-k-pivot_row_degree, max_deg) ;
2521 #endif /* NDEBUG */
2522
2523         /* === Add up set differences for each column ======================= */
2524
2525         DEBUG3 (("** Adding set differences phase. **\n")) ;
2526
2527         /* for each column in pivot row */
2528         rp = &A [pivot_row_start] ;
2529         rp_end = rp + pivot_row_length ;
2530         while (rp < rp_end)
2531         {
2532             /* get a column */
2533             col = *rp++ ;
2534             ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2535             hash = 0 ;
2536             cur_score = 0 ;
2537             cp = &A [Col [col].start] ;
2538             /* compact the column */
2539             new_cp = cp ;
2540             cp_end = cp + Col [col].length ;
2541
2542             DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ;
2543
2544             while (cp < cp_end)
2545             {
2546                 /* get a row */
2547                 row = *cp++ ;
2548                 ASSERT(row >= 0 && row < n_row) ;
2549                 row_mark = Row [row].shared2.mark ;
2550                 /* skip if dead */
2551                 if (ROW_IS_MARKED_DEAD (row_mark))
2552                 {
2553                     DEBUG4 ((" Row %d, dead\n", row)) ;
2554                     continue ;
2555                 }
2556                 DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark));
2557                 ASSERT (row_mark >= tag_mark) ;
2558                 /* compact the column */
2559                 *new_cp++ = row ;
2560                 /* compute hash function */
2561                 hash += row ;
2562                 /* add set difference */
2563                 cur_score += row_mark - tag_mark ;
2564                 /* integer overflow... */
2565                 cur_score = MIN (cur_score, n_col) ;
2566             }
2567
2568             /* recompute the column's length */
2569             Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
2570
2571             /* === Further mass elimination ================================= */
2572
2573             if (Col [col].length == 0)
2574             {
2575                 DEBUG4 (("further mass elimination. Col: %d\n", col)) ;
2576                 /* nothing left but the pivot row in this column */
2577                 KILL_PRINCIPAL_COL (col) ;
2578                 pivot_row_degree -= Col [col].shared1.thickness ;
2579                 ASSERT (pivot_row_degree >= 0) ;
2580                 /* order it */
2581                 Col [col].shared2.order = k ;
2582                 /* increment order count by column thickness */
2583                 k += Col [col].shared1.thickness ;
2584             }
2585             else
2586             {
2587                 /* === Prepare for supercolumn detection ==================== */
2588
2589                 DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ;
2590
2591                 /* save score so far */
2592                 Col [col].shared2.score = cur_score ;
2593
2594                 /* add column to hash table, for supercolumn detection */
2595                 hash %= n_col + 1 ;
2596
2597                 DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ;
2598                 ASSERT (((Int) hash) <= n_col) ;
2599
2600                 head_column = head [hash] ;
2601                 if (head_column > EMPTY)
2602                 {
2603                     /* degree list "hash" is non-empty, use prev (shared3) of */
2604                     /* first column in degree list as head of hash bucket */
2605                     first_col = Col [head_column].shared3.headhash ;
2606                     Col [head_column].shared3.headhash = col ;
2607                 }
2608                 else
2609                 {
2610                     /* degree list "hash" is empty, use head as hash bucket */
2611                     first_col = - (head_column + 2) ;
2612                     head [hash] = - (col + 2) ;
2613                 }
2614                 Col [col].shared4.hash_next = first_col ;
2615
2616                 /* save hash function in Col [col].shared3.hash */
2617                 Col [col].shared3.hash = (Int) hash ;
2618                 ASSERT (COL_IS_ALIVE (col)) ;
2619             }
2620         }
2621
2622         /* The approximate external column degree is now computed.  */
2623
2624         /* === Supercolumn detection ======================================== */
2625
2626         DEBUG3 (("** Supercolumn detection phase. **\n")) ;
2627
2628         detect_super_cols (
2629
2630 #ifndef NDEBUG
2631                 n_col, Row,
2632 #endif /* NDEBUG */
2633
2634                 Col, A, head, pivot_row_start, pivot_row_length) ;
2635
2636         /* === Kill the pivotal column ====================================== */
2637
2638         KILL_PRINCIPAL_COL (pivot_col) ;
2639
2640         /* === Clear mark =================================================== */
2641
2642         tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
2643
2644 #ifndef NDEBUG
2645         DEBUG3 (("check3\n")) ;
2646         debug_mark (n_row, Row, tag_mark, max_mark) ;
2647 #endif /* NDEBUG */
2648
2649         /* === Finalize the new pivot row, and column scores ================ */
2650
2651         DEBUG3 (("** Finalize scores phase. **\n")) ;
2652
2653         /* for each column in pivot row */
2654         rp = &A [pivot_row_start] ;
2655         /* compact the pivot row */
2656         new_rp = rp ;
2657         rp_end = rp + pivot_row_length ;
2658         while (rp < rp_end)
2659         {
2660             col = *rp++ ;
2661             /* skip dead columns */
2662             if (COL_IS_DEAD (col))
2663             {
2664                 continue ;
2665             }
2666             *new_rp++ = col ;
2667             /* add new pivot row to column */
2668             A [Col [col].start + (Col [col].length++)] = pivot_row ;
2669
2670             /* retrieve score so far and add on pivot row's degree. */
2671             /* (we wait until here for this in case the pivot */
2672             /* row's degree was reduced due to mass elimination). */
2673             cur_score = Col [col].shared2.score + pivot_row_degree ;
2674
2675             /* calculate the max possible score as the number of */
2676             /* external columns minus the 'k' value minus the */
2677             /* columns thickness */
2678             max_score = n_col - k - Col [col].shared1.thickness ;
2679
2680             /* make the score the external degree of the union-of-rows */
2681             cur_score -= Col [col].shared1.thickness ;
2682
2683             /* make sure score is less or equal than the max score */
2684             cur_score = MIN (cur_score, max_score) ;
2685             ASSERT (cur_score >= 0) ;
2686
2687             /* store updated score */
2688             Col [col].shared2.score = cur_score ;
2689
2690             /* === Place column back in degree list ========================= */
2691
2692             ASSERT (min_score >= 0) ;
2693             ASSERT (min_score <= n_col) ;
2694             ASSERT (cur_score >= 0) ;
2695             ASSERT (cur_score <= n_col) ;
2696             ASSERT (head [cur_score] >= EMPTY) ;
2697             next_col = head [cur_score] ;
2698             Col [col].shared4.degree_next = next_col ;
2699             Col [col].shared3.prev = EMPTY ;
2700             if (next_col != EMPTY)
2701             {
2702                 Col [next_col].shared3.prev = col ;
2703             }
2704             head [cur_score] = col ;
2705
2706             /* see if this score is less than current min */
2707             min_score = MIN (min_score, cur_score) ;
2708
2709         }
2710
2711 #ifndef NDEBUG
2712         debug_deg_lists (n_row, n_col, Row, Col, head,
2713                 min_score, n_col2-k, max_deg) ;
2714 #endif /* NDEBUG */
2715
2716         /* === Resurrect the new pivot row ================================== */
2717
2718         if (pivot_row_degree > 0)
2719         {
2720             /* update pivot row length to reflect any cols that were killed */
2721             /* during super-col detection and mass elimination */
2722             Row [pivot_row].start  = pivot_row_start ;
2723             Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
2724             ASSERT (Row [pivot_row].length > 0) ;
2725             Row [pivot_row].shared1.degree = pivot_row_degree ;
2726             Row [pivot_row].shared2.mark = 0 ;
2727             /* pivot row is no longer dead */
2728
2729             DEBUG1 (("Resurrect Pivot_row %d deg: %d\n",
2730                         pivot_row, pivot_row_degree)) ;
2731         }
2732     }
2733
2734     /* === All principal columns have now been ordered ====================== */
2735
2736     return (ngarbage) ;
2737 }
2738
2739
2740 /* ========================================================================== */
2741 /* === order_children ======================================================= */
2742 /* ========================================================================== */
2743
2744 /*
2745     The find_ordering routine has ordered all of the principal columns (the
2746     representatives of the supercolumns).  The non-principal columns have not
2747     yet been ordered.  This routine orders those columns by walking up the
2748     parent tree (a column is a child of the column which absorbed it).  The
2749     final permutation vector is then placed in p [0 ... n_col-1], with p [0]
2750     being the first column, and p [n_col-1] being the last.  It doesn't look
2751     like it at first glance, but be assured that this routine takes time linear
2752     in the number of columns.  Although not immediately obvious, the time
2753     taken by this routine is O (n_col), that is, linear in the number of
2754     columns.  Not user-callable.
2755 */
2756
2757 PRIVATE void order_children
2758 (
2759     /* === Parameters ======================================================= */
2760
2761     Int n_col,                  /* number of columns of A */
2762     Colamd_Col Col [],          /* of size n_col+1 */
2763     Int p []                    /* p [0 ... n_col-1] is the column permutation*/
2764 )
2765 {
2766     /* === Local variables ================================================== */
2767
2768     Int i ;                     /* loop counter for all columns */
2769     Int c ;                     /* column index */
2770     Int parent ;                /* index of column's parent */
2771     Int order ;                 /* column's order */
2772
2773     /* === Order each non-principal column ================================== */
2774
2775     for (i = 0 ; i < n_col ; i++)
2776     {
2777         /* find an un-ordered non-principal column */
2778         ASSERT (COL_IS_DEAD (i)) ;
2779         if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY)
2780         {
2781             parent = i ;
2782             /* once found, find its principal parent */
2783             do
2784             {
2785                 parent = Col [parent].shared1.parent ;
2786             } while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
2787
2788             /* now, order all un-ordered non-principal columns along path */
2789             /* to this parent.  collapse tree at the same time */
2790             c = i ;
2791             /* get order of parent */
2792             order = Col [parent].shared2.order ;
2793
2794             do
2795             {
2796                 ASSERT (Col [c].shared2.order == EMPTY) ;
2797
2798                 /* order this column */
2799                 Col [c].shared2.order = order++ ;
2800                 /* collaps tree */
2801                 Col [c].shared1.parent = parent ;
2802
2803                 /* get immediate parent of this column */
2804                 c = Col [c].shared1.parent ;
2805
2806                 /* continue until we hit an ordered column.  There are */
2807                 /* guarranteed not to be anymore unordered columns */
2808                 /* above an ordered column */
2809             } while (Col [c].shared2.order == EMPTY) ;
2810
2811             /* re-order the super_col parent to largest order for this group */
2812             Col [parent].shared2.order = order ;
2813         }
2814     }
2815
2816     /* === Generate the permutation ========================================= */
2817
2818     for (c = 0 ; c < n_col ; c++)
2819     {
2820         p [Col [c].shared2.order] = c ;
2821     }
2822 }
2823
2824
2825 /* ========================================================================== */
2826 /* === detect_super_cols ==================================================== */
2827 /* ========================================================================== */
2828
2829 /*
2830     Detects supercolumns by finding matches between columns in the hash buckets.
2831     Check amongst columns in the set A [row_start ... row_start + row_length-1].
2832     The columns under consideration are currently *not* in the degree lists,
2833     and have already been placed in the hash buckets.
2834
2835     The hash bucket for columns whose hash function is equal to h is stored
2836     as follows:
2837
2838         if head [h] is >= 0, then head [h] contains a degree list, so:
2839
2840                 head [h] is the first column in degree bucket h.
2841                 Col [head [h]].headhash gives the first column in hash bucket h.
2842
2843         otherwise, the degree list is empty, and:
2844
2845                 -(head [h] + 2) is the first column in hash bucket h.
2846
2847     For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
2848     column" pointer.  Col [c].shared3.hash is used instead as the hash number
2849     for that column.  The value of Col [c].shared4.hash_next is the next column
2850     in the same hash bucket.
2851
2852     Assuming no, or "few" hash collisions, the time taken by this routine is
2853     linear in the sum of the sizes (lengths) of each column whose score has
2854     just been computed in the approximate degree computation.
2855     Not user-callable.
2856 */
2857
2858 PRIVATE void detect_super_cols
2859 (
2860     /* === Parameters ======================================================= */
2861
2862 #ifndef NDEBUG
2863     /* these two parameters are only needed when debugging is enabled: */
2864     Int n_col,                  /* number of columns of A */
2865     Colamd_Row Row [],          /* of size n_row+1 */
2866 #endif /* NDEBUG */
2867
2868     Colamd_Col Col [],          /* of size n_col+1 */
2869     Int A [],                   /* row indices of A */
2870     Int head [],                /* head of degree lists and hash buckets */
2871     Int row_start,              /* pointer to set of columns to check */
2872     Int row_length              /* number of columns to check */
2873 )
2874 {
2875     /* === Local variables ================================================== */
2876
2877     Int hash ;                  /* hash value for a column */
2878     Int *rp ;                   /* pointer to a row */
2879     Int c ;                     /* a column index */
2880     Int super_c ;               /* column index of the column to absorb into */
2881     Int *cp1 ;                  /* column pointer for column super_c */
2882     Int *cp2 ;                  /* column pointer for column c */
2883     Int length ;                /* length of column super_c */
2884     Int prev_c ;                /* column preceding c in hash bucket */
2885     Int i ;                     /* loop counter */
2886     Int *rp_end ;               /* pointer to the end of the row */
2887     Int col ;                   /* a column index in the row to check */
2888     Int head_column ;           /* first column in hash bucket or degree list */
2889     Int first_col ;             /* first column in hash bucket */
2890
2891     /* === Consider each column in the row ================================== */
2892
2893     rp = &A [row_start] ;
2894     rp_end = rp + row_length ;
2895     while (rp < rp_end)
2896     {
2897         col = *rp++ ;
2898         if (COL_IS_DEAD (col))
2899         {
2900             continue ;
2901         }
2902
2903         /* get hash number for this column */
2904         hash = Col [col].shared3.hash ;
2905         ASSERT (hash <= n_col) ;
2906
2907         /* === Get the first column in this hash bucket ===================== */
2908
2909         head_column = head [hash] ;
2910         if (head_column > EMPTY)
2911         {
2912             first_col = Col [head_column].shared3.headhash ;
2913         }
2914         else
2915         {
2916             first_col = - (head_column + 2) ;
2917         }
2918
2919         /* === Consider each column in the hash bucket ====================== */
2920
2921         for (super_c = first_col ; super_c != EMPTY ;
2922             super_c = Col [super_c].shared4.hash_next)
2923         {
2924             ASSERT (COL_IS_ALIVE (super_c)) ;
2925             ASSERT (Col [super_c].shared3.hash == hash) ;
2926             length = Col [super_c].length ;
2927
2928             /* prev_c is the column preceding column c in the hash bucket */
2929             prev_c = super_c ;
2930
2931             /* === Compare super_c with all columns after it ================ */
2932
2933             for (c = Col [super_c].shared4.hash_next ;
2934                  c != EMPTY ; c = Col [c].shared4.hash_next)
2935             {
2936                 ASSERT (c != super_c) ;
2937                 ASSERT (COL_IS_ALIVE (c)) ;
2938                 ASSERT (Col [c].shared3.hash == hash) ;
2939
2940                 /* not identical if lengths or scores are different */
2941                 if (Col [c].length != length ||
2942                     Col [c].shared2.score != Col [super_c].shared2.score)
2943                 {
2944                     prev_c = c ;
2945                     continue ;
2946                 }
2947
2948                 /* compare the two columns */
2949                 cp1 = &A [Col [super_c].start] ;
2950                 cp2 = &A [Col [c].start] ;
2951
2952                 for (i = 0 ; i < length ; i++)
2953                 {
2954                     /* the columns are "clean" (no dead rows) */
2955                     ASSERT (ROW_IS_ALIVE (*cp1))  ;
2956                     ASSERT (ROW_IS_ALIVE (*cp2))  ;
2957                     /* row indices will same order for both supercols, */
2958                     /* no gather scatter nessasary */
2959                     if (*cp1++ != *cp2++)
2960                     {
2961                         break ;
2962                     }
2963                 }
2964
2965                 /* the two columns are different if the for-loop "broke" */
2966                 if (i != length)
2967                 {
2968                     prev_c = c ;
2969                     continue ;
2970                 }
2971
2972                 /* === Got it!  two columns are identical =================== */
2973
2974                 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
2975
2976                 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
2977                 Col [c].shared1.parent = super_c ;
2978                 KILL_NON_PRINCIPAL_COL (c) ;
2979                 /* order c later, in order_children() */
2980                 Col [c].shared2.order = EMPTY ;
2981                 /* remove c from hash bucket */
2982                 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
2983             }
2984         }
2985
2986         /* === Empty this hash bucket ======================================= */
2987
2988         if (head_column > EMPTY)
2989         {
2990             /* corresponding degree list "hash" is not empty */
2991             Col [head_column].shared3.headhash = EMPTY ;
2992         }
2993         else
2994         {
2995             /* corresponding degree list "hash" is empty */
2996             head [hash] = EMPTY ;
2997         }
2998     }
2999 }
3000
3001
3002 /* ========================================================================== */
3003 /* === garbage_collection =================================================== */
3004 /* ========================================================================== */
3005
3006 /*
3007     Defragments and compacts columns and rows in the workspace A.  Used when
3008     all avaliable memory has been used while performing row merging.  Returns
3009     the index of the first free position in A, after garbage collection.  The
3010     time taken by this routine is linear is the size of the array A, which is
3011     itself linear in the number of nonzeros in the input matrix.
3012     Not user-callable.
3013 */
3014
3015 PRIVATE Int garbage_collection  /* returns the new value of pfree */
3016 (
3017     /* === Parameters ======================================================= */
3018
3019     Int n_row,                  /* number of rows */
3020     Int n_col,                  /* number of columns */
3021     Colamd_Row Row [],          /* row info */
3022     Colamd_Col Col [],          /* column info */
3023     Int A [],                   /* A [0 ... Alen-1] holds the matrix */
3024     Int *pfree                  /* &A [0] ... pfree is in use */
3025 )
3026 {
3027     /* === Local variables ================================================== */
3028
3029     Int *psrc ;                 /* source pointer */
3030     Int *pdest ;                /* destination pointer */
3031     Int j ;                     /* counter */
3032     Int r ;                     /* a row index */
3033     Int c ;                     /* a column index */
3034     Int length ;                /* length of a row or column */
3035
3036 #ifndef NDEBUG
3037     Int debug_rows ;
3038     DEBUG2 (("Defrag..\n")) ;
3039     for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
3040     debug_rows = 0 ;
3041 #endif /* NDEBUG */
3042
3043     /* === Defragment the columns =========================================== */
3044
3045     pdest = &A[0] ;
3046     for (c = 0 ; c < n_col ; c++)
3047     {
3048         if (COL_IS_ALIVE (c))
3049         {
3050             psrc = &A [Col [c].start] ;
3051
3052             /* move and compact the column */
3053             ASSERT (pdest <= psrc) ;
3054             Col [c].start = (Int) (pdest - &A [0]) ;
3055             length = Col [c].length ;
3056             for (j = 0 ; j < length ; j++)
3057             {
3058                 r = *psrc++ ;
3059                 if (ROW_IS_ALIVE (r))
3060                 {
3061                     *pdest++ = r ;
3062                 }
3063             }
3064             Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
3065         }
3066     }
3067
3068     /* === Prepare to defragment the rows =================================== */
3069
3070     for (r = 0 ; r < n_row ; r++)
3071     {
3072         if (ROW_IS_DEAD (r) || (Row [r].length == 0))
3073         {
3074             /* This row is already dead, or is of zero length.  Cannot compact
3075              * a row of zero length, so kill it.  NOTE: in the current version,
3076              * there are no zero-length live rows.  Kill the row (for the first
3077              * time, or again) just to be safe. */
3078             KILL_ROW (r) ;
3079         }
3080         else
3081         {
3082             /* save first column index in Row [r].shared2.first_column */
3083             psrc = &A [Row [r].start] ;
3084             Row [r].shared2.first_column = *psrc ;
3085             ASSERT (ROW_IS_ALIVE (r)) ;
3086             /* flag the start of the row with the one's complement of row */
3087             *psrc = ONES_COMPLEMENT (r) ;
3088 #ifndef NDEBUG
3089             debug_rows++ ;
3090 #endif /* NDEBUG */
3091         }
3092     }
3093
3094     /* === Defragment the rows ============================================== */
3095
3096     psrc = pdest ;
3097     while (psrc < pfree)
3098     {
3099         /* find a negative number ... the start of a row */
3100         if (*psrc++ < 0)
3101         {
3102             psrc-- ;
3103             /* get the row index */
3104             r = ONES_COMPLEMENT (*psrc) ;
3105             ASSERT (r >= 0 && r < n_row) ;
3106             /* restore first column index */
3107             *psrc = Row [r].shared2.first_column ;
3108             ASSERT (ROW_IS_ALIVE (r)) ;
3109             ASSERT (Row [r].length > 0) ;
3110             /* move and compact the row */
3111             ASSERT (pdest <= psrc) ;
3112             Row [r].start = (Int) (pdest - &A [0]) ;
3113             length = Row [r].length ;
3114             for (j = 0 ; j < length ; j++)
3115             {
3116                 c = *psrc++ ;
3117                 if (COL_IS_ALIVE (c))
3118                 {
3119                     *pdest++ = c ;
3120                 }
3121             }
3122             Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
3123             ASSERT (Row [r].length > 0) ;
3124 #ifndef NDEBUG
3125             debug_rows-- ;
3126 #endif /* NDEBUG */
3127         }
3128     }
3129     /* ensure we found all the rows */
3130     ASSERT (debug_rows == 0) ;
3131
3132     /* === Return the new value of pfree ==================================== */
3133
3134     return ((Int) (pdest - &A [0])) ;
3135 }
3136
3137
3138 /* ========================================================================== */
3139 /* === clear_mark =========================================================== */
3140 /* ========================================================================== */
3141
3142 /*
3143     Clears the Row [].shared2.mark array, and returns the new tag_mark.
3144     Return value is the new tag_mark.  Not user-callable.
3145 */
3146
3147 PRIVATE Int clear_mark  /* return the new value for tag_mark */
3148 (
3149     /* === Parameters ======================================================= */
3150
3151     Int tag_mark,       /* new value of tag_mark */
3152     Int max_mark,       /* max allowed value of tag_mark */
3153
3154     Int n_row,          /* number of rows in A */
3155     Colamd_Row Row []   /* Row [0 ... n_row-1].shared2.mark is set to zero */
3156 )
3157 {
3158     /* === Local variables ================================================== */
3159
3160     Int r ;
3161
3162     if (tag_mark <= 0 || tag_mark >= max_mark)
3163     {
3164         for (r = 0 ; r < n_row ; r++)
3165         {
3166             if (ROW_IS_ALIVE (r))
3167             {
3168                 Row [r].shared2.mark = 0 ;
3169             }
3170         }
3171         tag_mark = 1 ;
3172     }
3173
3174     return (tag_mark) ;
3175 }
3176
3177
3178 /* ========================================================================== */
3179 /* === print_report ========================================================= */
3180 /* ========================================================================== */
3181
3182 PRIVATE void print_report
3183 (
3184     char *method,
3185     Int stats [COLAMD_STATS]
3186 )
3187 {
3188
3189     Int i1, i2, i3 ;
3190
3191     PRINTF (("\n%s version %d.%d, %s: ", method,
3192             COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ;
3193
3194     if (!stats)
3195     {
3196         PRINTF (("No statistics available.\n")) ;
3197         return ;
3198     }
3199
3200     i1 = stats [COLAMD_INFO1] ;
3201     i2 = stats [COLAMD_INFO2] ;
3202     i3 = stats [COLAMD_INFO3] ;
3203
3204     if (stats [COLAMD_STATUS] >= 0)
3205     {
3206         PRINTF (("OK.  ")) ;
3207     }
3208     else
3209     {
3210         PRINTF (("ERROR.  ")) ;
3211     }
3212
3213     switch (stats [COLAMD_STATUS])
3214     {
3215
3216         case COLAMD_OK_BUT_JUMBLED:
3217
3218             PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
3219
3220             PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n",
3221             method, i3)) ;
3222
3223             PRINTF(("%s: last seen duplicate or out-of-order row index:   %d\n",
3224             method, INDEX (i2))) ;
3225
3226             PRINTF(("%s: last seen in column:                             %d",
3227             method, INDEX (i1))) ;
3228
3229             /* no break - fall through to next case instead */
3230
3231         case COLAMD_OK:
3232
3233             PRINTF(("\n")) ;
3234
3235             PRINTF(("%s: number of dense or empty rows ignored:           %d\n",
3236             method, stats [COLAMD_DENSE_ROW])) ;
3237
3238             PRINTF(("%s: number of dense or empty columns ignored:        %d\n",
3239             method, stats [COLAMD_DENSE_COL])) ;
3240
3241             PRINTF(("%s: number of garbage collections performed:         %d\n",
3242             method, stats [COLAMD_DEFRAG_COUNT])) ;
3243             break ;
3244
3245         case COLAMD_ERROR_A_not_present:
3246
3247             PRINTF(("Array A (row indices of matrix) not present.\n")) ;
3248             break ;
3249
3250         case COLAMD_ERROR_p_not_present:
3251
3252             PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
3253             break ;
3254
3255         case COLAMD_ERROR_nrow_negative:
3256
3257             PRINTF(("Invalid number of rows (%d).\n", i1)) ;
3258             break ;
3259
3260         case COLAMD_ERROR_ncol_negative:
3261
3262             PRINTF(("Invalid number of columns (%d).\n", i1)) ;
3263             break ;
3264
3265         case COLAMD_ERROR_nnz_negative:
3266
3267             PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ;
3268             break ;
3269
3270         case COLAMD_ERROR_p0_nonzero:
3271
3272             PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1));
3273             break ;
3274
3275         case COLAMD_ERROR_A_too_small:
3276
3277             PRINTF(("Array A too small.\n")) ;
3278             PRINTF(("        Need Alen >= %d, but given only Alen = %d.\n",
3279             i1, i2)) ;
3280             break ;
3281
3282         case COLAMD_ERROR_col_length_negative:
3283
3284             PRINTF
3285             (("Column %d has a negative number of nonzero entries (%d).\n",
3286             INDEX (i1), i2)) ;
3287             break ;
3288
3289         case COLAMD_ERROR_row_index_out_of_bounds:
3290
3291             PRINTF
3292             (("Row index (row %d) out of bounds (%d to %d) in column %d.\n",
3293             INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ;
3294             break ;
3295
3296         case COLAMD_ERROR_out_of_memory:
3297
3298             PRINTF(("Out of memory.\n")) ;
3299             break ;
3300
3301         /* v2.4: internal-error case deleted */
3302     }
3303 }
3304
3305
3306
3307
3308 /* ========================================================================== */
3309 /* === colamd debugging routines ============================================ */
3310 /* ========================================================================== */
3311
3312 /* When debugging is disabled, the remainder of this file is ignored. */
3313
3314 #ifndef NDEBUG
3315
3316
3317 /* ========================================================================== */
3318 /* === debug_structures ===================================================== */
3319 /* ========================================================================== */
3320
3321 /*
3322     At this point, all empty rows and columns are dead.  All live columns
3323     are "clean" (containing no dead rows) and simplicial (no supercolumns
3324     yet).  Rows may contain dead columns, but all live rows contain at
3325     least one live column.
3326 */
3327
3328 PRIVATE void debug_structures
3329 (
3330     /* === Parameters ======================================================= */
3331
3332     Int n_row,
3333     Int n_col,
3334     Colamd_Row Row [],
3335     Colamd_Col Col [],
3336     Int A [],
3337     Int n_col2
3338 )
3339 {
3340     /* === Local variables ================================================== */
3341
3342     Int i ;
3343     Int c ;
3344     Int *cp ;
3345     Int *cp_end ;
3346     Int len ;
3347     Int score ;
3348     Int r ;
3349     Int *rp ;
3350     Int *rp_end ;
3351     Int deg ;
3352
3353     /* === Check A, Row, and Col ============================================ */
3354
3355     for (c = 0 ; c < n_col ; c++)
3356     {
3357         if (COL_IS_ALIVE (c))
3358         {
3359             len = Col [c].length ;
3360             score = Col [c].shared2.score ;
3361             DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
3362             ASSERT (len > 0) ;
3363             ASSERT (score >= 0) ;
3364             ASSERT (Col [c].shared1.thickness == 1) ;
3365             cp = &A [Col [c].start] ;
3366             cp_end = cp + len ;
3367             while (cp < cp_end)
3368             {
3369                 r = *cp++ ;
3370                 ASSERT (ROW_IS_ALIVE (r)) ;
3371             }
3372         }
3373         else
3374         {
3375             i = Col [c].shared2.order ;
3376             ASSERT (i >= n_col2 && i < n_col) ;
3377         }
3378     }
3379
3380     for (r = 0 ; r < n_row ; r++)
3381     {
3382         if (ROW_IS_ALIVE (r))
3383         {
3384             i = 0 ;
3385             len = Row [r].length ;
3386             deg = Row [r].shared1.degree ;
3387             ASSERT (len > 0) ;
3388             ASSERT (deg > 0) ;
3389             rp = &A [Row [r].start] ;
3390             rp_end = rp + len ;
3391             while (rp < rp_end)
3392             {
3393                 c = *rp++ ;
3394                 if (COL_IS_ALIVE (c))
3395                 {
3396                     i++ ;
3397                 }
3398             }
3399             ASSERT (i > 0) ;
3400         }
3401     }
3402 }
3403