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. NumLoopsLoopExtractorPass24 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 668 #if 0 /* by mao */ 669 #include <limits.h> 670 #include <math.h> 671 672 #ifdef MATLAB_MEX_FILE 673 #include "mex.h" 674 #include "matrix.h" 675 #endif /* MATLAB_MEX_FILE */ 676 677 #if !defined (NPRINT) || !defined (NDEBUG) 678 #include <stdio.h> 679 #endif 680 681 #ifndef NULL 682 #define NULL ((void *) 0) 683 #endif 684 #endif 685 686 /* ========================================================================== */ 687 /* === int or UF_long ======================================================= */ 688 /* ========================================================================== */ 689 690 #if 0 /* by mao */ 691 /* define UF_long */ 692 #include "UFconfig.h" 693 #endif 694 695 #ifdef DLONG 696 697 #define Int UF_long 698 #define ID UF_long_id 699 #define Int_MAX UF_long_max 700 701 #define COLAMD_recommended colamd_l_recommended 702 #define COLAMD_set_defaults colamd_l_set_defaults 703 #define COLAMD_MAIN colamd_l 704 #define SYMAMD_MAIN symamd_l 705 #define COLAMD_report colamd_l_report 706 #define SYMAMD_report symamd_l_report 707 708 #else 709 710 #define Int int 711 #define ID "%d" 712 #define Int_MAX INT_MAX 713 714 #define COLAMD_recommended colamd_recommended 715 #define COLAMD_set_defaults colamd_set_defaults 716 #define COLAMD_MAIN colamd 717 #define SYMAMD_MAIN symamd 718 #define COLAMD_report colamd_report 719 #define SYMAMD_report symamd_report 720 721 #endif 722 723 /* ========================================================================== */ 724 /* === Row and Column structures ============================================ */ 725 /* ========================================================================== */ 726 727 /* User code that makes use of the colamd/symamd routines need not directly */ 728 /* reference these structures. They are used only for colamd_recommended. */ 729 730 typedef struct Colamd_Col_struct 731 { 732 Int start ; /* index for A of first row in this column, or DEAD */ 733 /* if column is dead */ 734 Int length ; /* number of rows in this column */ 735 union 736 { 737 Int thickness ; /* number of original columns represented by this */ 738 /* col, if the column is alive */ 739 Int parent ; /* parent in parent tree super-column structure, if */ 740 /* the column is dead */ 741 } shared1 ; 742 union 743 { 744 Int score ; /* the score used to maintain heap, if col is alive */ 745 Int order ; /* pivot ordering of this column, if col is dead */ 746 } shared2 ; 747 union 748 { 749 Int headhash ; /* head of a hash bucket, if col is at the head of */ 750 /* a degree list */ 751 Int hash ; /* hash value, if col is not in a degree list */ 752 Int prev ; /* previous column in degree list, if col is in a */ 753 /* degree list (but not at the head of a degree list) */ 754 } shared3 ; 755 union 756 { 757 Int degree_next ; /* next column, if col is in a degree list */ 758 Int hash_next ; /* next column, if col is in a hash list */ 759 } shared4 ; 760 761 } Colamd_Col ; 762 763 typedef struct Colamd_Row_struct 764 { 765 Int start ; /* index for A of first col in this row */ 766 Int length ; /* number of principal columns in this row */ 767 union 768 { 769 Int degree ; /* number of principal & non-principal columns in row */ 770 Int p ; /* used as a row pointer in init_rows_cols () */ 771 } shared1 ; 772 union 773 { 774 Int mark ; /* for computing set differences and marking dead rows*/ 775 Int first_column ;/* first column in row (used in garbage collection) */ 776 } shared2 ; 777 778 } Colamd_Row ; 779 780 /* ========================================================================== */ 781 /* === Definitions ========================================================== */ 782 /* ========================================================================== */ 783 784 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */ 785 #define PUBLIC 786 #define PRIVATE static 787 788 #define DENSE_DEGREE(alpha,n) \ 789 ((Int) MAX (16.0, (alpha) * sqrt ((double) (n)))) 790 791 #define MAX(a,b) (((a) > (b)) ? (a) : (b)) 792 #define MIN(a,b) (((a) < (b)) ? (a) : (b)) 793 794 #define ONES_COMPLEMENT(r) (-(r)-1) 795 796 /* -------------------------------------------------------------------------- */ 797 /* Change for version 2.1: define TRUE and FALSE only if not yet defined */ 798 /* -------------------------------------------------------------------------- */ 799 800 #ifndef TRUE 801 #define TRUE (1) 802 #endif 803 804 #ifndef FALSE 805 #define FALSE (0) 806 #endif 807 808 /* -------------------------------------------------------------------------- */ 809 810 #define EMPTY (-1) 811 812 /* Row and column status */ 813 #define ALIVE (0) 814 #define DEAD (-1) 815 816 /* Column status */ 817 #define DEAD_PRINCIPAL (-1) 818 #define DEAD_NON_PRINCIPAL (-2) 819 820 /* Macros for row and column status update and checking. */ 821 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) 822 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) 823 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) 824 #define COL_IS_DEAD(c) (Col [c].start < ALIVE) 825 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) 826 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) 827 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } 828 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } 829 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } 830 831 /* ========================================================================== */ 832 /* === Colamd reporting mechanism =========================================== */ 833 /* ========================================================================== */ 834 835 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS) 836 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */ 837 #define INDEX(i) ((i)+1) 838 #else 839 /* In C, matrices are 0-based and indices are reported as such in *_report */ 840 #define INDEX(i) (i) 841 #endif 842 843 /* All output goes through the PRINTF macro. */ 844 #define PRINTF(params) { if (colamd_printf != NULL) (void) colamd_printf params ; } 845 846 /* ========================================================================== */ 847 /* === Prototypes of PRIVATE routines ======================================= */ 848 /* ========================================================================== */ 849 850 PRIVATE Int init_rows_cols 851 ( 852 Int n_row, 853 Int n_col, 854 Colamd_Row Row [], 855 Colamd_Col Col [], 856 Int A [], 857 Int p [], 858 Int stats [COLAMD_STATS] 859 ) ; 860 861 PRIVATE void init_scoring 862 ( 863 Int n_row, 864 Int n_col, 865 Colamd_Row Row [], 866 Colamd_Col Col [], 867 Int A [], 868 Int head [], 869 double knobs [COLAMD_KNOBS], 870 Int *p_n_row2, 871 Int *p_n_col2, 872 Int *p_max_deg 873 ) ; 874 875 PRIVATE Int find_ordering 876 ( 877 Int n_row, 878 Int n_col, 879 Int Alen, 880 Colamd_Row Row [], 881 Colamd_Col Col [], 882 Int A [], 883 Int head [], 884 Int n_col2, 885 Int max_deg, 886 Int pfree, 887 Int aggressive 888 ) ; 889 890 PRIVATE void order_children 891 ( 892 Int n_col, 893 Colamd_Col Col [], 894 Int p [] 895 ) ; 896 897 PRIVATE void detect_super_cols 898 ( 899 900 #ifndef NDEBUG 901 Int n_col, 902 Colamd_Row Row [], 903 #endif /* NDEBUG */ 904 905 Colamd_Col Col [], 906 Int A [], 907 Int head [], 908 Int row_start, 909 Int row_length 910 ) ; 911 912 PRIVATE Int garbage_collection 913 ( 914 Int n_row, 915 Int n_col, 916 Colamd_Row Row [], 917 Colamd_Col Col [], 918 Int A [], 919 Int *pfree 920 ) ; 921 922 PRIVATE Int clear_mark 923 ( 924 Int tag_mark, 925 Int max_mark, 926 Int n_row, 927 Colamd_Row Row [] 928 ) ; 929 930 PRIVATE void print_report 931 ( 932 char *method, 933 Int stats [COLAMD_STATS] 934 ) ; 935 936 /* ========================================================================== */ 937 /* === Debugging prototypes and definitions ================================= */ 938 /* ========================================================================== */ 939 940 #ifndef NDEBUG 941 942 #if 0 /* by mao */ 943 #include <assert.h> 944 #endif 945 946 /* colamd_debug is the *ONLY* global variable, and is only */ 947 /* present when debugging */ 948 949 PRIVATE Int colamd_debug = 0 ; /* debug print level */ 950 951 #define DEBUG0(params) { PRINTF (params) ; } 952 #define DEBUG1(params) { if (colamd_debug >= 1) PRINTF (params) ; } 953 #define DEBUG2(params) { if (colamd_debug >= 2) PRINTF (params) ; } 954 #define DEBUG3(params) { if (colamd_debug >= 3) PRINTF (params) ; } 955 #define DEBUG4(params) { if (colamd_debug >= 4) PRINTF (params) ; } 956 957 #if 0 /* by mao */ 958 #ifdef MATLAB_MEX_FILE 959 #define ASSERT(expression) (mxAssert ((expression), "")) 960 #else 961 #define ASSERT(expression) (assert (expression)) 962 #endif /* MATLAB_MEX_FILE */ 963 #else 964 #define ASSERT xassert 965 #endif 966 967 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */ 968 ( 969 char *method 970 ) ; 971 972 PRIVATE void debug_deg_lists 973 ( 974 Int n_row, 975 Int n_col, 976 Colamd_Row Row [], 977 Colamd_Col Col [], 978 Int head [], 979 Int min_score, 980 Int should, 981 Int max_deg 982 ) ; 983 984 PRIVATE void debug_mark 985 ( 986 Int n_row, 987 Colamd_Row Row [], 988 Int tag_mark, 989 Int max_mark 990 ) ; 991 992 PRIVATE void debug_matrix 993 ( 994 Int n_row, 995 Int n_col, 996 Colamd_Row Row [], 997 Colamd_Col Col [], 998 Int A [] 999 ) ; 1000 1001 PRIVATE void debug_structures 1002 ( 1003 Int n_row, 1004 Int n_col, 1005 Colamd_Row Row [], 1006 Colamd_Col Col [], 1007 Int A [], 1008 Int n_col2 1009 ) ; 1010 1011 #else /* NDEBUG */ 1012 1013 /* === No debugging ========================================================= */ 1014 1015 #define DEBUG0(params) ; 1016 #define DEBUG1(params) ; 1017 #define DEBUG2(params) ; 1018 #define DEBUG3(params) ; 1019 #define DEBUG4(params) ; 1020 1021 #define ASSERT(expression) 1022 1023 #endif /* NDEBUG */ 1024 1025 /* ========================================================================== */ 1026 /* === USER-CALLABLE ROUTINES: ============================================== */ 1027 /* ========================================================================== */ 1028 1029 /* ========================================================================== */ 1030 /* === colamd_recommended =================================================== */ 1031 /* ========================================================================== */ 1032 1033 /* 1034 The colamd_recommended routine returns the suggested size for Alen. This 1035 value has been determined to provide good balance between the number of 1036 garbage collections and the memory requirements for colamd. If any 1037 argument is negative, or if integer overflow occurs, a 0 is returned as an 1038 error condition. 2*nnz space is required for the row and column 1039 indices of the matrix. COLAMD_C (n_col) + COLAMD_R (n_row) space is 1040 required for the Col and Row arrays, respectively, which are internal to 1041 colamd (roughly 6*n_col + 4*n_row). An additional n_col space is the 1042 minimal amount of "elbow room", and nnz/5 more space is recommended for 1043 run time efficiency. 1044 1045 Alen is approximately 2.2*nnz + 7*n_col + 4*n_row + 10. 1046 1047 This function is not needed when using symamd. 1048 */ 1049 1050 /* add two values of type size_t, and check for integer overflow */ 1051 static size_t t_add (size_t a, size_t b, int *ok) 1052 { 1053 (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ; 1054 return ((*ok) ? (a + b) : 0) ; 1055 } 1056 1057 /* compute a*k where k is a small integer, and check for integer overflow */ 1058 static size_t t_mult (size_t a, size_t k, int *ok) 1059 { 1060 size_t i, s = 0 ; 1061 for (i = 0 ; i < k ; i++) 1062 { 1063 s = t_add (s, a, ok) ; 1064 } 1065 return (s) ; 1066 } 1067 1068 /* size of the Col and Row structures */ 1069 #define COLAMD_C(n_col,ok) \ 1070 ((t_mult (t_add (n_col, 1, ok), sizeof (Colamd_Col), ok) / sizeof (Int))) 1071 1072 #define COLAMD_R(n_row,ok) \ 1073 ((t_mult (t_add (n_row, 1, ok), sizeof (Colamd_Row), ok) / sizeof (Int))) 1074 1075 1076 PUBLIC size_t COLAMD_recommended /* returns recommended value of Alen. */ 1077 ( 1078 /* === Parameters ======================================================= */ 1079 1080 Int nnz, /* number of nonzeros in A */ 1081 Int n_row, /* number of rows in A */ 1082 Int n_col /* number of columns in A */ 1083 ) 1084 { 1085 size_t s, c, r ; 1086 int ok = TRUE ; 1087 if (nnz < 0 || n_row < 0 || n_col < 0) 1088 { 1089 return (0) ; 1090 } 1091 s = t_mult (nnz, 2, &ok) ; /* 2*nnz */ 1092 c = COLAMD_C (n_col, &ok) ; /* size of column structures */ 1093 r = COLAMD_R (n_row, &ok) ; /* size of row structures */ 1094 s = t_add (s, c, &ok) ; 1095 s = t_add (s, r, &ok) ; 1096 s = t_add (s, n_col, &ok) ; /* elbow room */ 1097 s = t_add (s, nnz/5, &ok) ; /* elbow room */ 1098 ok = ok && (s < Int_MAX) ; 1099 return (ok ? s : 0) ; 1100 } 1101 1102 1103 /* ========================================================================== */ 1104 /* === colamd_set_defaults ================================================== */ 1105 /* ========================================================================== */ 1106 1107 /* 1108 The colamd_set_defaults routine sets the default values of the user- 1109 controllable parameters for colamd and symamd: 1110 1111 Colamd: rows with more than max (16, knobs [0] * sqrt (n_col)) 1112 entries are removed prior to ordering. Columns with more than 1113 max (16, knobs [1] * sqrt (MIN (n_row,n_col))) entries are removed 1114 prior to ordering, and placed last in the output column ordering. 1115 1116 Symamd: Rows and columns with more than max (16, knobs [0] * sqrt (n)) 1117 entries are removed prior to ordering, and placed last in the 1118 output ordering. 1119 1120 knobs [0] dense row control 1121 1122 knobs [1] dense column control 1123 1124 knobs [2] if nonzero, do aggresive absorption 1125 1126 knobs [3..19] unused, but future versions might use this 1127 1128 */ 1129 1130 PUBLIC void COLAMD_set_defaults 1131 ( 1132 /* === Parameters ======================================================= */ 1133 1134 double knobs [COLAMD_KNOBS] /* knob array */ 1135 ) 1136 { 1137 /* === Local variables ================================================== */ 1138 1139 Int i ; 1140 1141 if (!knobs) 1142 { 1143 return ; /* no knobs to initialize */ 1144 } 1145 for (i = 0 ; i < COLAMD_KNOBS ; i++) 1146 { 1147 knobs [i] = 0 ; 1148 } 1149 knobs [COLAMD_DENSE_ROW] = 10 ; 1150 knobs [COLAMD_DENSE_COL] = 10 ; 1151 knobs [COLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/ 1152 } 1153 1154 1155 /* ========================================================================== */ 1156 /* === symamd =============================================================== */ 1157 /* ========================================================================== */ 1158 1159 PUBLIC Int SYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */ 1160 ( 1161 /* === Parameters ======================================================= */ 1162 1163 Int n, /* number of rows and columns of A */ 1164 Int A [], /* row indices of A */ 1165 Int p [], /* column pointers of A */ 1166 Int perm [], /* output permutation, size n+1 */ 1167 double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */ 1168 Int stats [COLAMD_STATS], /* output statistics and error codes */ 1169 void * (*allocate) (size_t, size_t), 1170 /* pointer to calloc (ANSI C) or */ 1171 /* mxCalloc (for MATLAB mexFunction) */ 1172 void (*release) (void *) 1173 /* pointer to free (ANSI C) or */ 1174 /* mxFree (for MATLAB mexFunction) */ 1175 ) 1176 { 1177 /* === Local variables ================================================== */ 1178 1179 Int *count ; /* length of each column of M, and col pointer*/ 1180 Int *mark ; /* mark array for finding duplicate entries */ 1181 Int *M ; /* row indices of matrix M */ 1182 size_t Mlen ; /* length of M */ 1183 Int n_row ; /* number of rows in M */ 1184 Int nnz ; /* number of entries in A */ 1185 Int i ; /* row index of A */ 1186 Int j ; /* column index of A */ 1187 Int k ; /* row index of M */ 1188 Int mnz ; /* number of nonzeros in M */ 1189 Int pp ; /* index into a column of A */ 1190 Int last_row ; /* last row seen in the current column */ 1191 Int length ; /* number of nonzeros in a column */ 1192 1193 double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */ 1194 double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */ 1195 1196 #ifndef NDEBUG 1197 colamd_get_debug ("symamd") ; 1198 #endif /* NDEBUG */ 1199 1200 /* === Check the input arguments ======================================== */ 1201 1202 if (!stats) 1203 { 1204 DEBUG0 (("symamd: stats not present\n")) ; 1205 return (FALSE) ; 1206 } 1207 for (i = 0 ; i < COLAMD_STATS ; i++) 1208 { 1209 stats [i] = 0 ; 1210 } 1211 stats [COLAMD_STATUS] = COLAMD_OK ; 1212 stats [COLAMD_INFO1] = -1 ; 1213 stats [COLAMD_INFO2] = -1 ; 1214 1215 if (!A) 1216 { 1217 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; 1218 DEBUG0 (("symamd: A not present\n")) ; 1219 return (FALSE) ; 1220 } 1221 1222 if (!p) /* p is not present */ 1223 { 1224 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; 1225 DEBUG0 (("symamd: p not present\n")) ; 1226 return (FALSE) ; 1227 } 1228 1229 if (n < 0) /* n must be >= 0 */ 1230 { 1231 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; 1232 stats [COLAMD_INFO1] = n ; 1233 DEBUG0 (("symamd: n negative %d\n", n)) ; 1234 return (FALSE) ; 1235 } 1236 1237 nnz = p [n] ; 1238 if (nnz < 0) /* nnz must be >= 0 */ 1239 { 1240 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; 1241 stats [COLAMD_INFO1] = nnz ; 1242 DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ; 1243 return (FALSE) ; 1244 } 1245 1246 if (p [0] != 0) 1247 { 1248 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; 1249 stats [COLAMD_INFO1] = p [0] ; 1250 DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ; 1251 return (FALSE) ; 1252 } 1253 1254 /* === If no knobs, set default knobs =================================== */ 1255 1256 if (!knobs) 1257 { 1258 COLAMD_set_defaults (default_knobs) ; 1259 knobs = default_knobs ; 1260 } 1261 1262 /* === Allocate count and mark ========================================== */ 1263 1264 count = (Int *) ((*allocate) (n+1, sizeof (Int))) ; 1265 if (!count) 1266 { 1267 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; 1268 DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ; 1269 return (FALSE) ; 1270 } 1271 1272 mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ; 1273 if (!mark) 1274 { 1275 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; 1276 (*release) ((void *) count) ; 1277 DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ; 1278 return (FALSE) ; 1279 } 1280 1281 /* === Compute column counts of M, check if A is valid ================== */ 1282 1283 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ 1284 1285 for (i = 0 ; i < n ; i++) 1286 { 1287 mark [i] = -1 ; 1288 } 1289 1290 for (j = 0 ; j < n ; j++) 1291 { 1292 last_row = -1 ; 1293 1294 length = p [j+1] - p [j] ; 1295 if (length < 0) 1296 { 1297 /* column pointers must be non-decreasing */ 1298 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; 1299 stats [COLAMD_INFO1] = j ; 1300 stats [COLAMD_INFO2] = length ; 1301 (*release) ((void *) count) ; 1302 (*release) ((void *) mark) ; 1303 DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ; 1304 return (FALSE) ; 1305 } 1306 1307 for (pp = p [j] ; pp < p [j+1] ; pp++) 1308 { 1309 i = A [pp] ; 1310 if (i < 0 || i >= n) 1311 { 1312 /* row index i, in column j, is out of bounds */ 1313 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; 1314 stats [COLAMD_INFO1] = j ; 1315 stats [COLAMD_INFO2] = i ; 1316 stats [COLAMD_INFO3] = n ; 1317 (*release) ((void *) count) ; 1318 (*release) ((void *) mark) ; 1319 DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ; 1320 return (FALSE) ; 1321 } 1322 1323 if (i <= last_row || mark [i] == j) 1324 { 1325 /* row index is unsorted or repeated (or both), thus col */ 1326 /* is jumbled. This is a notice, not an error condition. */ 1327 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; 1328 stats [COLAMD_INFO1] = j ; 1329 stats [COLAMD_INFO2] = i ; 1330 (stats [COLAMD_INFO3]) ++ ; 1331 DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ; 1332 } 1333 1334 if (i > j && mark [i] != j) 1335 { 1336 /* row k of M will contain column indices i and j */ 1337 count [i]++ ; 1338 count [j]++ ; 1339 } 1340 1341 /* mark the row as having been seen in this column */ 1342 mark [i] = j ; 1343 1344 last_row = i ; 1345 } 1346 } 1347 1348 /* v2.4: removed free(mark) */ 1349 1350 /* === Compute column pointers of M ===================================== */ 1351 1352 /* use output permutation, perm, for column pointers of M */ 1353 perm [0] = 0 ; 1354 for (j = 1 ; j <= n ; j++) 1355 { 1356 perm [j] = perm [j-1] + count [j-1] ; 1357 } 1358 for (j = 0 ; j < n ; j++) 1359 { 1360 count [j] = perm [j] ; 1361 } 1362 1363 /* === Construct M ====================================================== */ 1364 1365 mnz = perm [n] ; 1366 n_row = mnz / 2 ; 1367 Mlen = COLAMD_recommended (mnz, n_row, n) ; 1368 M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ; 1369 DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %g\n", 1370 n_row, n, mnz, (double) Mlen)) ; 1371 1372 if (!M) 1373 { 1374 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; 1375 (*release) ((void *) count) ; 1376 (*release) ((void *) mark) ; 1377 DEBUG0 (("symamd: allocate M (size %g) failed\n", (double) Mlen)) ; 1378 return (FALSE) ; 1379 } 1380 1381 k = 0 ; 1382 1383 if (stats [COLAMD_STATUS] == COLAMD_OK) 1384 { 1385 /* Matrix is OK */ 1386 for (j = 0 ; j < n ; j++) 1387 { 1388 ASSERT (p [j+1] - p [j] >= 0) ; 1389 for (pp = p [j] ; pp < p [j+1] ; pp++) 1390 { 1391 i = A [pp] ; 1392 ASSERT (i >= 0 && i < n) ; 1393 if (i > j) 1394 { 1395 /* row k of M contains column indices i and j */ 1396 M [count [i]++] = k ; 1397 M [count [j]++] = k ; 1398 k++ ; 1399 } 1400 } 1401 } 1402 } 1403 else 1404 { 1405 /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */ 1406 DEBUG0 (("symamd: Duplicates in A.\n")) ; 1407 for (i = 0 ; i < n ; i++) 1408 { 1409 mark [i] = -1 ; 1410 } 1411 for (j = 0 ; j < n ; j++) 1412 { 1413 ASSERT (p [j+1] - p [j] >= 0) ; 1414 for (pp = p [j] ; pp < p [j+1] ; pp++) 1415 { 1416 i = A [pp] ; 1417 ASSERT (i >= 0 && i < n) ; 1418 if (i > j && mark [i] != j) 1419 { 1420 /* row k of M contains column indices i and j */ 1421 M [count [i]++] = k ; 1422 M [count [j]++] = k ; 1423 k++ ; 1424 mark [i] = j ; 1425 } 1426 } 1427 } 1428 /* v2.4: free(mark) moved below */ 1429 } 1430 1431 /* count and mark no longer needed */ 1432 (*release) ((void *) count) ; 1433 (*release) ((void *) mark) ; /* v2.4: free (mark) moved here */ 1434 ASSERT (k == n_row) ; 1435 1436 /* === Adjust the knobs for M =========================================== */ 1437 1438 for (i = 0 ; i < COLAMD_KNOBS ; i++) 1439 { 1440 cknobs [i] = knobs [i] ; 1441 } 1442 1443 /* there are no dense rows in M */ 1444 cknobs [COLAMD_DENSE_ROW] = -1 ; 1445 cknobs [COLAMD_DENSE_COL] = knobs [COLAMD_DENSE_ROW] ; 1446 1447 /* === Order the columns of M =========================================== */ 1448 1449 /* v2.4: colamd cannot fail here, so the error check is removed */ 1450 (void) COLAMD_MAIN (n_row, n, (Int) Mlen, M, perm, cknobs, stats) ; 1451 1452 /* Note that the output permutation is now in perm */ 1453 1454 /* === get the statistics for symamd from colamd ======================== */ 1455 1456 /* a dense column in colamd means a dense row and col in symamd */ 1457 stats [COLAMD_DENSE_ROW] = stats [COLAMD_DENSE_COL] ; 1458 1459 /* === Free M =========================================================== */ 1460 1461 (*release) ((void *) M) ; 1462 DEBUG0 (("symamd: done.\n")) ; 1463 return (TRUE) ; 1464 1465 } 1466 1467 /* ========================================================================== */ 1468 /* === colamd =============================================================== */ 1469 /* ========================================================================== */ 1470 1471 /* 1472 The colamd routine computes a column ordering Q of a sparse matrix 1473 A such that the LU factorization P(AQ) = LU remains sparse, where P is 1474 selected via partial pivoting. The routine can also be viewed as 1475 providing a permutation Q such that the Cholesky factorization 1476 (AQ)'(AQ) = LL' remains sparse. 1477 */ 1478 1479 PUBLIC Int COLAMD_MAIN /* returns TRUE if successful, FALSE otherwise*/ 1480 ( 1481 /* === Parameters ======================================================= */ 1482 1483 Int n_row, /* number of rows in A */ 1484 Int n_col, /* number of columns in A */ 1485 Int Alen, /* length of A */ 1486 Int A [], /* row indices of A */ 1487 Int p [], /* pointers to columns in A */ 1488 double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */ 1489 Int stats [COLAMD_STATS] /* output statistics and error codes */ 1490 ) 1491 { 1492 /* === Local variables ================================================== */ 1493 1494 Int i ; /* loop index */ 1495 Int nnz ; /* nonzeros in A */ 1496 size_t Row_size ; /* size of Row [], in integers */ 1497 size_t Col_size ; /* size of Col [], in integers */ 1498 size_t need ; /* minimum required length of A */ 1499 Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ 1500 Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ 1501 Int n_col2 ; /* number of non-dense, non-empty columns */ 1502 Int n_row2 ; /* number of non-dense, non-empty rows */ 1503 Int ngarbage ; /* number of garbage collections performed */ 1504 Int max_deg ; /* maximum row degree */ 1505 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ 1506 Int aggressive ; /* do aggressive absorption */ 1507 int ok ; 1508 1509 #ifndef NDEBUG 1510 colamd_get_debug ("colamd") ; 1511 #endif /* NDEBUG */ 1512 1513 /* === Check the input arguments ======================================== */ 1514 1515 if (!stats) 1516 { 1517 DEBUG0 (("colamd: stats not present\n")) ; 1518 return (FALSE) ; 1519 } 1520 for (i = 0 ; i < COLAMD_STATS ; i++) 1521 { 1522 stats [i] = 0 ; 1523 } 1524 stats [COLAMD_STATUS] = COLAMD_OK ; 1525 stats [COLAMD_INFO1] = -1 ; 1526 stats [COLAMD_INFO2] = -1 ; 1527 1528 if (!A) /* A is not present */ 1529 { 1530 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; 1531 DEBUG0 (("colamd: A not present\n")) ; 1532 return (FALSE) ; 1533 } 1534 1535 if (!p) /* p is not present */ 1536 { 1537 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; 1538 DEBUG0 (("colamd: p not present\n")) ; 1539 return (FALSE) ; 1540 } 1541 1542 if (n_row < 0) /* n_row must be >= 0 */ 1543 { 1544 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; 1545 stats [COLAMD_INFO1] = n_row ; 1546 DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; 1547 return (FALSE) ; 1548 } 1549 1550 if (n_col < 0) /* n_col must be >= 0 */ 1551 { 1552 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; 1553 stats [COLAMD_INFO1] = n_col ; 1554 DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; 1555 return (FALSE) ; 1556 } 1557 1558 nnz = p [n_col] ; 1559 if (nnz < 0) /* nnz must be >= 0 */ 1560 { 1561 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; 1562 stats [COLAMD_INFO1] = nnz ; 1563 DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; 1564 return (FALSE) ; 1565 } 1566 1567 if (p [0] != 0) 1568 { 1569 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; 1570 stats [COLAMD_INFO1] = p [0] ; 1571 DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; 1572 return (FALSE) ; 1573 } 1574 1575 /* === If no knobs, set default knobs =================================== */ 1576 1577 if (!knobs) 1578 { 1579 COLAMD_set_defaults (default_knobs) ; 1580 knobs = default_knobs ; 1581 } 1582 1583 aggressive = (knobs [COLAMD_AGGRESSIVE] != FALSE) ; 1584 1585 /* === Allocate the Row and Col arrays from array A ===================== */ 1586 1587 ok = TRUE ; 1588 Col_size = COLAMD_C (n_col, &ok) ; /* size of Col array of structs */ 1589 Row_size = COLAMD_R (n_row, &ok) ; /* size of Row array of structs */ 1590 1591 /* need = 2*nnz + n_col + Col_size + Row_size ; */ 1592 need = t_mult (nnz, 2, &ok) ; 1593 need = t_add (need, n_col, &ok) ; 1594 need = t_add (need, Col_size, &ok) ; 1595 need = t_add (need, Row_size, &ok) ; 1596 1597 if (!ok || need > (size_t) Alen || need > Int_MAX) 1598 { 1599 /* not enough space in array A to perform the ordering */ 1600 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; 1601 stats [COLAMD_INFO1] = need ; 1602 stats [COLAMD_INFO2] = Alen ; 1603 DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); 1604 return (FALSE) ; 1605 } 1606 1607 Alen -= Col_size + Row_size ; 1608 Col = (Colamd_Col *) &A [Alen] ; 1609 Row = (Colamd_Row *) &A [Alen + Col_size] ; 1610 1611 /* === Construct the row and column data structures ===================== */ 1612 1613 if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) 1614 { 1615 /* input matrix is invalid */ 1616 DEBUG0 (("colamd: Matrix invalid\n")) ; 1617 return (FALSE) ; 1618 } 1619 1620 /* === Initialize scores, kill dense rows/columns ======================= */ 1621 1622 init_scoring (n_row, n_col, Row, Col, A, p, knobs, 1623 &n_row2, &n_col2, &max_deg) ; 1624 1625 /* === Order the supercolumns =========================================== */ 1626 1627 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, 1628 n_col2, max_deg, 2*nnz, aggressive) ; 1629 1630 /* === Order the non-principal columns ================================== */ 1631 1632 order_children (n_col, Col, p) ; 1633 1634 /* === Return statistics in stats ======================================= */ 1635 1636 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; 1637 stats [COLAMD_DENSE_COL] = n_col - n_col2 ; 1638 stats [COLAMD_DEFRAG_COUNT] = ngarbage ; 1639 DEBUG0 (("colamd: done.\n")) ; 1640 return (TRUE) ; 1641 } 1642 1643 1644 /* ========================================================================== */ 1645 /* === colamd_report ======================================================== */ 1646 /* ========================================================================== */ 1647 1648 PUBLIC void COLAMD_report 1649 ( 1650 Int stats [COLAMD_STATS] 1651 ) 1652 { 1653 print_report ("colamd", stats) ; 1654 } 1655 1656 1657 /* ========================================================================== */ 1658 /* === symamd_report ======================================================== */ 1659 /* ========================================================================== */ 1660 1661 PUBLIC void SYMAMD_report 1662 ( 1663 Int stats [COLAMD_STATS] 1664 ) 1665 { 1666 print_report ("symamd", stats) ; 1667 } 1668 1669 1670 1671 /* ========================================================================== */ 1672 /* === NON-USER-CALLABLE ROUTINES: ========================================== */ 1673 /* ========================================================================== */ 1674 1675 /* There are no user-callable routines beyond this point in the file */ 1676 1677 1678 /* ========================================================================== */ 1679 /* === init_rows_cols ======================================================= */ 1680 /* ========================================================================== */ 1681 1682 /* 1683 Takes the column form of the matrix in A and creates the row form of the 1684 matrix. Also, row and column attributes are stored in the Col and Row 1685 structs. If the columns are un-sorted or contain duplicate row indices, 1686 this routine will also sort and remove duplicate row indices from the 1687 column form of the matrix. Returns FALSE if the matrix is invalid, 1688 TRUE otherwise. Not user-callable. 1689 */ 1690 1691 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */ 1692 ( 1693 /* === Parameters ======================================================= */ 1694 1695 Int n_row, /* number of rows of A */ 1696 Int n_col, /* number of columns of A */ 1697 Colamd_Row Row [], /* of size n_row+1 */ 1698 Colamd_Col Col [], /* of size n_col+1 */ 1699 Int A [], /* row indices of A, of size Alen */ 1700 Int p [], /* pointers to columns in A, of size n_col+1 */ 1701 Int stats [COLAMD_STATS] /* colamd statistics */ 1702 ) 1703 { 1704 /* === Local variables ================================================== */ 1705 1706 Int col ; /* a column index */ 1707 Int row ; /* a row index */ 1708 Int *cp ; /* a column pointer */ 1709 Int *cp_end ; /* a pointer to the end of a column */ 1710 Int *rp ; /* a row pointer */ 1711 Int *rp_end ; /* a pointer to the end of a row */ 1712 Int last_row ; /* previous row */ 1713 1714 /* === Initialize columns, and check column pointers ==================== */ 1715 1716 for (col = 0 ; col < n_col ; col++) 1717 { 1718 Col [col].start = p [col] ; 1719 Col [col].length = p [col+1] - p [col] ; 1720 1721 if (Col [col].length < 0) 1722 { 1723 /* column pointers must be non-decreasing */ 1724 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; 1725 stats [COLAMD_INFO1] = col ; 1726 stats [COLAMD_INFO2] = Col [col].length ; 1727 DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; 1728 return (FALSE) ; 1729 } 1730 1731 Col [col].shared1.thickness = 1 ; 1732 Col [col].shared2.score = 0 ; 1733 Col [col].shared3.prev = EMPTY ; 1734 Col [col].shared4.degree_next = EMPTY ; 1735 } 1736 1737 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ 1738 1739 /* === Scan columns, compute row degrees, and check row indices ========= */ 1740 1741 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ 1742 1743 for (row = 0 ; row < n_row ; row++) 1744 { 1745 Row [row].length = 0 ; 1746 Row [row].shared2.mark = -1 ; 1747 } 1748 1749 for (col = 0 ; col < n_col ; col++) 1750 { 1751 last_row = -1 ; 1752 1753 cp = &A [p [col]] ; 1754 cp_end = &A [p [col+1]] ; 1755 1756 while (cp < cp_end) 1757 { 1758 row = *cp++ ; 1759 1760 /* make sure row indices within range */ 1761 if (row < 0 || row >= n_row) 1762 { 1763 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; 1764 stats [COLAMD_INFO1] = col ; 1765 stats [COLAMD_INFO2] = row ; 1766 stats [COLAMD_INFO3] = n_row ; 1767 DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; 1768 return (FALSE) ; 1769 } 1770 1771 if (row <= last_row || Row [row].shared2.mark == col) 1772 { 1773 /* row index are unsorted or repeated (or both), thus col */ 1774 /* is jumbled. This is a notice, not an error condition. */ 1775 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; 1776 stats [COLAMD_INFO1] = col ; 1777 stats [COLAMD_INFO2] = row ; 1778 (stats [COLAMD_INFO3]) ++ ; 1779 DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); 1780 } 1781 1782 if (Row [row].shared2.mark != col) 1783 { 1784 Row [row].length++ ; 1785 } 1786 else 1787 { 1788 /* this is a repeated entry in the column, */ 1789 /* it will be removed */ 1790 Col [col].length-- ; 1791 } 1792 1793 /* mark the row as having been seen in this column */ 1794 Row [row].shared2.mark = col ; 1795 1796 last_row = row ; 1797 } 1798 } 1799 1800 /* === Compute row pointers ============================================= */ 1801 1802 /* row form of the matrix starts directly after the column */ 1803 /* form of matrix in A */ 1804 Row [0].start = p [n_col] ; 1805 Row [0].shared1.p = Row [0].start ; 1806 Row [0].shared2.mark = -1 ; 1807 for (row = 1 ; row < n_row ; row++) 1808 { 1809 Row [row].start = Row [row-1].start + Row [row-1].length ; 1810 Row [row].shared1.p = Row [row].start ; 1811 Row [row].shared2.mark = -1 ; 1812 } 1813 1814 /* === Create row form ================================================== */ 1815 1816 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 1817 { 1818 /* if cols jumbled, watch for repeated row indices */ 1819 for (col = 0 ; col < n_col ; col++) 1820 { 1821 cp = &A [p [col]] ; 1822 cp_end = &A [p [col+1]] ; 1823 while (cp < cp_end) 1824 { 1825 row = *cp++ ; 1826 if (Row [row].shared2.mark != col) 1827 { 1828 A [(Row [row].shared1.p)++] = col ; 1829 Row [row].shared2.mark = col ; 1830 } 1831 } 1832 } 1833 } 1834 else 1835 { 1836 /* if cols not jumbled, we don't need the mark (this is faster) */ 1837 for (col = 0 ; col < n_col ; col++) 1838 { 1839 cp = &A [p [col]] ; 1840 cp_end = &A [p [col+1]] ; 1841 while (cp < cp_end) 1842 { 1843 A [(Row [*cp++].shared1.p)++] = col ; 1844 } 1845 } 1846 } 1847 1848 /* === Clear the row marks and set row degrees ========================== */ 1849 1850 for (row = 0 ; row < n_row ; row++) 1851 { 1852 Row [row].shared2.mark = 0 ; 1853 Row [row].shared1.degree = Row [row].length ; 1854 } 1855 1856 /* === See if we need to re-create columns ============================== */ 1857 1858 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 1859 { 1860 DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; 1861 1862 #ifndef NDEBUG 1863 /* make sure column lengths are correct */ 1864 for (col = 0 ; col < n_col ; col++) 1865 { 1866 p [col] = Col [col].length ; 1867 } 1868 for (row = 0 ; row < n_row ; row++) 1869 { 1870 rp = &A [Row [row].start] ; 1871 rp_end = rp + Row [row].length ; 1872 while (rp < rp_end) 1873 { 1874 p [*rp++]-- ; 1875 } 1876 } 1877 for (col = 0 ; col < n_col ; col++) 1878 { 1879 ASSERT (p [col] == 0) ; 1880 } 1881 /* now p is all zero (different than when debugging is turned off) */ 1882 #endif /* NDEBUG */ 1883 1884 /* === Compute col pointers ========================================= */ 1885 1886 /* col form of the matrix starts at A [0]. */ 1887 /* Note, we may have a gap between the col form and the row */ 1888 /* form if there were duplicate entries, if so, it will be */ 1889 /* removed upon the first garbage collection */ 1890 Col [0].start = 0 ; 1891 p [0] = Col [0].start ; 1892 for (col = 1 ; col < n_col ; col++) 1893 { 1894 /* note that the lengths here are for pruned columns, i.e. */ 1895 /* no duplicate row indices will exist for these columns */ 1896 Col [col].start = Col [col-1].start + Col [col-1].length ; 1897 p [col] = Col [col].start ; 1898 } 1899 1900 /* === Re-create col form =========================================== */ 1901 1902 for (row = 0 ; row < n_row ; row++) 1903 { 1904 rp = &A [Row [row].start] ; 1905 rp_end = rp + Row [row].length ; 1906 while (rp < rp_end) 1907 { 1908 A [(p [*rp++])++] = row ; 1909 } 1910 } 1911 } 1912 1913 /* === Done. Matrix is not (or no longer) jumbled ====================== */ 1914 1915 return (TRUE) ; 1916 } 1917 1918 1919 /* ========================================================================== */ 1920 /* === init_scoring ========================================================= */ 1921 /* ========================================================================== */ 1922 1923 /* 1924 Kills dense or empty columns and rows, calculates an initial score for 1925 each column, and places all columns in the degree lists. Not user-callable. 1926 */ 1927 1928 PRIVATE void init_scoring 1929 ( 1930 /* === Parameters ======================================================= */ 1931 1932 Int n_row, /* number of rows of A */ 1933 Int n_col, /* number of columns of A */ 1934 Colamd_Row Row [], /* of size n_row+1 */ 1935 Colamd_Col Col [], /* of size n_col+1 */ 1936 Int A [], /* column form and row form of A */ 1937 Int head [], /* of size n_col+1 */ 1938 double knobs [COLAMD_KNOBS],/* parameters */ 1939 Int *p_n_row2, /* number of non-dense, non-empty rows */ 1940 Int *p_n_col2, /* number of non-dense, non-empty columns */ 1941 Int *p_max_deg /* maximum row degree */ 1942 ) 1943 { 1944 /* === Local variables ================================================== */ 1945 1946 Int c ; /* a column index */ 1947 Int r, row ; /* a row index */ 1948 Int *cp ; /* a column pointer */ 1949 Int deg ; /* degree of a row or column */ 1950 Int *cp_end ; /* a pointer to the end of a column */ 1951 Int *new_cp ; /* new column pointer */ 1952 Int col_length ; /* length of pruned column */ 1953 Int score ; /* current column score */ 1954 Int n_col2 ; /* number of non-dense, non-empty columns */ 1955 Int n_row2 ; /* number of non-dense, non-empty rows */ 1956 Int dense_row_count ; /* remove rows with more entries than this */ 1957 Int dense_col_count ; /* remove cols with more entries than this */ 1958 Int min_score ; /* smallest column score */ 1959 Int max_deg ; /* maximum row degree */ 1960 Int next_col ; /* Used to add to degree list.*/ 1961 1962 #ifndef NDEBUG 1963 Int debug_count ; /* debug only. */ 1964 #endif /* NDEBUG */ 1965 1966 /* === Extract knobs ==================================================== */ 1967 1968 /* Note: if knobs contains a NaN, this is undefined: */ 1969 if (knobs [COLAMD_DENSE_ROW] < 0) 1970 { 1971 /* only remove completely dense rows */ 1972 dense_row_count = n_col-1 ; 1973 } 1974 else 1975 { 1976 dense_row_count = DENSE_DEGREE (knobs [COLAMD_DENSE_ROW], n_col) ; 1977 } 1978 if (knobs [COLAMD_DENSE_COL] < 0) 1979 { 1980 /* only remove completely dense columns */ 1981 dense_col_count = n_row-1 ; 1982 } 1983 else 1984 { 1985 dense_col_count = 1986 DENSE_DEGREE (knobs [COLAMD_DENSE_COL], MIN (n_row, n_col)) ; 1987 } 1988 1989 DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; 1990 max_deg = 0 ; 1991 n_col2 = n_col ; 1992 n_row2 = n_row ; 1993 1994 /* === Kill empty columns =============================================== */ 1995 1996 /* Put the empty columns at the end in their natural order, so that LU */ 1997 /* factorization can proceed as far as possible. */ 1998 for (c = n_col-1 ; c >= 0 ; c--) 1999 { 2000 deg = Col [c].length ; 2001 if (deg == 0) 2002 { 2003 /* this is a empty column, kill and order it last */ 2004 Col [c].shared2.order = --n_col2 ; 2005 KILL_PRINCIPAL_COL (c) ; 2006 } 2007 } 2008 DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; 2009 2010 /* === Kill dense columns =============================================== */ 2011 2012 /* Put the dense columns at the end, in their natural order */ 2013 for (c = n_col-1 ; c >= 0 ; c--) 2014 { 2015 /* skip any dead columns */ 2016 if (COL_IS_DEAD (c)) 2017 { 2018 continue ; 2019 } 2020 deg = Col [c].length ; 2021 if (deg > dense_col_count) 2022 { 2023 /* this is a dense column, kill and order it last */ 2024 Col [c].shared2.order = --n_col2 ; 2025 /* decrement the row degrees */ 2026 cp = &A [Col [c].start] ; 2027 cp_end = cp + Col [c].length ; 2028 while (cp < cp_end) 2029 { 2030 Row [*cp++].shared1.degree-- ; 2031 } 2032 KILL_PRINCIPAL_COL (c) ; 2033 } 2034 } 2035 DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; 2036 2037 /* === Kill dense and empty rows ======================================== */ 2038 2039 for (r = 0 ; r < n_row ; r++) 2040 { 2041 deg = Row [r].shared1.degree ; 2042 ASSERT (deg >= 0 && deg <= n_col) ; 2043 if (deg > dense_row_count || deg == 0) 2044 { 2045 /* kill a dense or empty row */ 2046 KILL_ROW (r) ; 2047 --n_row2 ; 2048 } 2049 else 2050 { 2051 /* keep track of max degree of remaining rows */ 2052 max_deg = MAX (max_deg, deg) ; 2053 } 2054 } 2055 DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; 2056 2057 /* === Compute initial column scores ==================================== */ 2058 2059 /* At this point the row degrees are accurate. They reflect the number */ 2060 /* of "live" (non-dense) columns in each row. No empty rows exist. */ 2061 /* Some "live" columns may contain only dead rows, however. These are */ 2062 /* pruned in the code below. */ 2063 2064 /* now find the initial matlab score for each column */ 2065 for (c = n_col-1 ; c >= 0 ; c--) 2066 { 2067 /* skip dead column */ 2068 if (COL_IS_DEAD (c)) 2069 { 2070 continue ; 2071 } 2072 score = 0 ; 2073 cp = &A [Col [c].start] ; 2074 new_cp = cp ; 2075 cp_end = cp + Col [c].length ; 2076 while (cp < cp_end) 2077 { 2078 /* get a row */ 2079 row = *cp++ ; 2080 /* skip if dead */ 2081 if (ROW_IS_DEAD (row)) 2082 { 2083 continue ; 2084 } 2085 /* compact the column */ 2086 *new_cp++ = row ; 2087 /* add row's external degree */ 2088 score += Row [row].shared1.degree - 1 ; 2089 /* guard against integer overflow */ 2090 score = MIN (score, n_col) ; 2091 } 2092 /* determine pruned column length */ 2093 col_length = (Int) (new_cp - &A [Col [c].start]) ; 2094 if (col_length == 0) 2095 { 2096 /* a newly-made null column (all rows in this col are "dense" */ 2097 /* and have already been killed) */ 2098 DEBUG2 (("Newly null killed: %d\n", c)) ; 2099 Col [c].shared2.order = --n_col2 ; 2100 KILL_PRINCIPAL_COL (c) ; 2101 } 2102 else 2103 { 2104 /* set column length and set score */ 2105 ASSERT (score >= 0) ; 2106 ASSERT (score <= n_col) ; 2107 Col [c].length = col_length ; 2108 Col [c].shared2.score = score ; 2109 } 2110 } 2111 DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", 2112 n_col-n_col2)) ; 2113 2114 /* At this point, all empty rows and columns are dead. All live columns */ 2115 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ 2116 /* yet). Rows may contain dead columns, but all live rows contain at */ 2117 /* least one live column. */ 2118 2119 #ifndef NDEBUG 2120 debug_structures (n_row, n_col, Row, Col, A, n_col2) ; 2121 #endif /* NDEBUG */ 2122 2123 /* === Initialize degree lists ========================================== */ 2124 2125 #ifndef NDEBUG 2126 debug_count = 0 ; 2127 #endif /* NDEBUG */ 2128 2129 /* clear the hash buckets */ 2130 for (c = 0 ; c <= n_col ; c++) 2131 { 2132 head [c] = EMPTY ; 2133 } 2134 min_score = n_col ; 2135 /* place in reverse order, so low column indices are at the front */ 2136 /* of the lists. This is to encourage natural tie-breaking */ 2137 for (c = n_col-1 ; c >= 0 ; c--) 2138 { 2139 /* only add principal columns to degree lists */ 2140 if (COL_IS_ALIVE (c)) 2141 { 2142 DEBUG4 (("place %d score %d minscore %d ncol %d\n", 2143 c, Col [c].shared2.score, min_score, n_col)) ; 2144 2145 /* === Add columns score to DList =============================== */ 2146 2147 score = Col [c].shared2.score ; 2148 2149 ASSERT (min_score >= 0) ; 2150 ASSERT (min_score <= n_col) ; 2151 ASSERT (score >= 0) ; 2152 ASSERT (score <= n_col) ; 2153 ASSERT (head [score] >= EMPTY) ; 2154 2155 /* now add this column to dList at proper score location */ 2156 next_col = head [score] ; 2157 Col [c].shared3.prev = EMPTY ; 2158 Col [c].shared4.degree_next = next_col ; 2159 2160 /* if there already was a column with the same score, set its */ 2161 /* previous pointer to this new column */ 2162 if (next_col != EMPTY) 2163 { 2164 Col [next_col].shared3.prev = c ; 2165 } 2166 head [score] = c ; 2167 2168 /* see if this score is less than current min */ 2169 min_score = MIN (min_score, score) ; 2170 2171 #ifndef NDEBUG 2172 debug_count++ ; 2173 #endif /* NDEBUG */ 2174 2175 } 2176 } 2177 2178 #ifndef NDEBUG 2179 DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n", 2180 debug_count, n_col, n_col-debug_count)) ; 2181 ASSERT (debug_count == n_col2) ; 2182 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; 2183 #endif /* NDEBUG */ 2184 2185 /* === Return number of remaining columns, and max row degree =========== */ 2186 2187 *p_n_col2 = n_col2 ; 2188 *p_n_row2 = n_row2 ; 2189 *p_max_deg = max_deg ; 2190 } 2191 2192 2193 /* ========================================================================== */ 2194 /* === find_ordering ======================================================== */ 2195 /* ========================================================================== */ 2196 2197 /* 2198 Order the principal columns of the supercolumn form of the matrix 2199 (no supercolumns on input). Uses a minimum approximate column minimum 2200 degree ordering method. Not user-callable. 2201 */ 2202 2203 PRIVATE Int find_ordering /* return the number of garbage collections */ 2204 ( 2205 /* === Parameters ======================================================= */ 2206 2207 Int n_row, /* number of rows of A */ 2208 Int n_col, /* number of columns of A */ 2209 Int Alen, /* size of A, 2*nnz + n_col or larger */ 2210 Colamd_Row Row [], /* of size n_row+1 */ 2211 Colamd_Col Col [], /* of size n_col+1 */ 2212 Int A [], /* column form and row form of A */ 2213 Int head [], /* of size n_col+1 */ 2214 Int n_col2, /* Remaining columns to order */ 2215 Int max_deg, /* Maximum row degree */ 2216 Int pfree, /* index of first free slot (2*nnz on entry) */ 2217 Int aggressive 2218 ) 2219 { 2220 /* === Local variables ================================================== */ 2221 2222 Int k ; /* current pivot ordering step */ 2223 Int pivot_col ; /* current pivot column */ 2224 Int *cp ; /* a column pointer */ 2225 Int *rp ; /* a row pointer */ 2226 Int pivot_row ; /* current pivot row */ 2227 Int *new_cp ; /* modified column pointer */ 2228 Int *new_rp ; /* modified row pointer */ 2229 Int pivot_row_start ; /* pointer to start of pivot row */ 2230 Int pivot_row_degree ; /* number of columns in pivot row */ 2231 Int pivot_row_length ; /* number of supercolumns in pivot row */ 2232 Int pivot_col_score ; /* score of pivot column */ 2233 Int needed_memory ; /* free space needed for pivot row */ 2234 Int *cp_end ; /* pointer to the end of a column */ 2235 Int *rp_end ; /* pointer to the end of a row */ 2236 Int row ; /* a row index */ 2237 Int col ; /* a column index */ 2238 Int max_score ; /* maximum possible score */ 2239 Int cur_score ; /* score of current column */ 2240 unsigned Int hash ; /* hash value for supernode detection */ 2241 Int head_column ; /* head of hash bucket */ 2242 Int first_col ; /* first column in hash bucket */ 2243 Int tag_mark ; /* marker value for mark array */ 2244 Int row_mark ; /* Row [row].shared2.mark */ 2245 Int set_difference ; /* set difference size of row with pivot row */ 2246 Int min_score ; /* smallest column score */ 2247 Int col_thickness ; /* "thickness" (no. of columns in a supercol) */ 2248 Int max_mark ; /* maximum value of tag_mark */ 2249 Int pivot_col_thickness ; /* number of columns represented by pivot col */ 2250 Int prev_col ; /* Used by Dlist operations. */ 2251 Int next_col ; /* Used by Dlist operations. */ 2252 Int ngarbage ; /* number of garbage collections performed */ 2253 2254 #ifndef NDEBUG 2255 Int debug_d ; /* debug loop counter */ 2256 Int debug_step = 0 ; /* debug loop counter */ 2257 #endif /* NDEBUG */ 2258 2259 /* === Initialization and clear mark ==================================== */ 2260 2261 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */ 2262 tag_mark = clear_mark (0, max_mark, n_row, Row) ; 2263 min_score = 0 ; 2264 ngarbage = 0 ; 2265 DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; 2266 2267 /* === Order the columns ================================================ */ 2268 2269 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) 2270 { 2271 2272 #ifndef NDEBUG 2273 if (debug_step % 100 == 0) 2274 { 2275 DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; 2276 } 2277 else 2278 { 2279 DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; 2280 } 2281 debug_step++ ; 2282 debug_deg_lists (n_row, n_col, Row, Col, head, 2283 min_score, n_col2-k, max_deg) ; 2284 debug_matrix (n_row, n_col, Row, Col, A) ; 2285 #endif /* NDEBUG */ 2286 2287 /* === Select pivot column, and order it ============================ */ 2288 2289 /* make sure degree list isn't empty */ 2290 ASSERT (min_score >= 0) ; 2291 ASSERT (min_score <= n_col) ; 2292 ASSERT (head [min_score] >= EMPTY) ; 2293 2294 #ifndef NDEBUG 2295 for (debug_d = 0 ; debug_d < min_score ; debug_d++) 2296 { 2297 ASSERT (head [debug_d] == EMPTY) ; 2298 } 2299 #endif /* NDEBUG */ 2300 2301 /* get pivot column from head of minimum degree list */ 2302 while (head [min_score] == EMPTY && min_score < n_col) 2303 { 2304 min_score++ ; 2305 } 2306 pivot_col = head [min_score] ; 2307 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; 2308 next_col = Col [pivot_col].shared4.degree_next ; 2309 head [min_score] = next_col ; 2310 if (next_col != EMPTY) 2311 { 2312 Col [next_col].shared3.prev = EMPTY ; 2313 } 2314 2315 ASSERT (COL_IS_ALIVE (pivot_col)) ; 2316 2317 /* remember score for defrag check */ 2318 pivot_col_score = Col [pivot_col].shared2.score ; 2319 2320 /* the pivot column is the kth column in the pivot order */ 2321 Col [pivot_col].shared2.order = k ; 2322 2323 /* increment order count by column thickness */ 2324 pivot_col_thickness = Col [pivot_col].shared1.thickness ; 2325 k += pivot_col_thickness ; 2326 ASSERT (pivot_col_thickness > 0) ; 2327 DEBUG3 (("Pivot col: %d thick %d\n", pivot_col, pivot_col_thickness)) ; 2328 2329 /* === Garbage_collection, if necessary ============================= */ 2330 2331 needed_memory = MIN (pivot_col_score, n_col - k) ; 2332 if (pfree + needed_memory >= Alen) 2333 { 2334 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; 2335 ngarbage++ ; 2336 /* after garbage collection we will have enough */ 2337 ASSERT (pfree + needed_memory < Alen) ; 2338 /* garbage collection has wiped out the Row[].shared2.mark array */ 2339 tag_mark = clear_mark (0, max_mark, n_row, Row) ; 2340 2341 #ifndef NDEBUG 2342 debug_matrix (n_row, n_col, Row, Col, A) ; 2343 #endif /* NDEBUG */ 2344 } 2345 2346 /* === Compute pivot row pattern ==================================== */ 2347 2348 /* get starting location for this new merged row */ 2349 pivot_row_start = pfree ; 2350 2351 /* initialize new row counts to zero */ 2352 pivot_row_degree = 0 ; 2353 2354 /* tag pivot column as having been visited so it isn't included */ 2355 /* in merged pivot row */ 2356 Col [pivot_col].shared1.thickness = -pivot_col_thickness ; 2357 2358 /* pivot row is the union of all rows in the pivot column pattern */ 2359 cp = &A [Col [pivot_col].start] ; 2360 cp_end = cp + Col [pivot_col].length ; 2361 while (cp < cp_end) 2362 { 2363 /* get a row */ 2364 row = *cp++ ; 2365 DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; 2366 /* skip if row is dead */ 2367 if (ROW_IS_ALIVE (row)) 2368 { 2369 rp = &A [Row [row].start] ; 2370 rp_end = rp + Row [row].length ; 2371 while (rp < rp_end) 2372 { 2373 /* get a column */ 2374 col = *rp++ ; 2375 /* add the column, if alive and untagged */ 2376 col_thickness = Col [col].shared1.thickness ; 2377 if (col_thickness > 0 && COL_IS_ALIVE (col)) 2378 { 2379 /* tag column in pivot row */ 2380 Col [col].shared1.thickness = -col_thickness ; 2381 ASSERT (pfree < Alen) ; 2382 /* place column in pivot row */ 2383 A [pfree++] = col ; 2384 pivot_row_degree += col_thickness ; 2385 } 2386 } 2387 } 2388 } 2389 2390 /* clear tag on pivot column */ 2391 Col [pivot_col].shared1.thickness = pivot_col_thickness ; 2392 max_deg = MAX (max_deg, pivot_row_degree) ; 2393 2394 #ifndef NDEBUG 2395 DEBUG3 (("check2\n")) ; 2396 debug_mark (n_row, Row, tag_mark, max_mark) ; 2397 #endif /* NDEBUG */ 2398 2399 /* === Kill all rows used to construct pivot row ==================== */ 2400 2401 /* also kill pivot row, temporarily */ 2402 cp = &A [Col [pivot_col].start] ; 2403 cp_end = cp + Col [pivot_col].length ; 2404 while (cp < cp_end) 2405 { 2406 /* may be killing an already dead row */ 2407 row = *cp++ ; 2408 DEBUG3 (("Kill row in pivot col: %d\n", row)) ; 2409 KILL_ROW (row) ; 2410 } 2411 2412 /* === Select a row index to use as the new pivot row =============== */ 2413 2414 pivot_row_length = pfree - pivot_row_start ; 2415 if (pivot_row_length > 0) 2416 { 2417 /* pick the "pivot" row arbitrarily (first row in col) */ 2418 pivot_row = A [Col [pivot_col].start] ; 2419 DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; 2420 } 2421 else 2422 { 2423 /* there is no pivot row, since it is of zero length */ 2424 pivot_row = EMPTY ; 2425 ASSERT (pivot_row_length == 0) ; 2426 } 2427 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; 2428 2429 /* === Approximate degree computation =============================== */ 2430 2431 /* Here begins the computation of the approximate degree. The column */ 2432 /* score is the sum of the pivot row "length", plus the size of the */ 2433 /* set differences of each row in the column minus the pattern of the */ 2434 /* pivot row itself. The column ("thickness") itself is also */ 2435 /* excluded from the column score (we thus use an approximate */ 2436 /* external degree). */ 2437 2438 /* The time taken by the following code (compute set differences, and */ 2439 /* add them up) is proportional to the size of the data structure */ 2440 /* being scanned - that is, the sum of the sizes of each column in */ 2441 /* the pivot row. Thus, the amortized time to compute a column score */ 2442 /* is proportional to the size of that column (where size, in this */ 2443 /* context, is the column "length", or the number of row indices */ 2444 /* in that column). The number of row indices in a column is */ 2445 /* monotonically non-decreasing, from the length of the original */ 2446 /* column on input to colamd. */ 2447 2448 /* === Compute set differences ====================================== */ 2449 2450 DEBUG3 (("** Computing set differences phase. **\n")) ; 2451 2452 /* pivot row is currently dead - it will be revived later. */ 2453 2454 DEBUG3 (("Pivot row: ")) ; 2455 /* for each column in pivot row */ 2456 rp = &A [pivot_row_start] ; 2457 rp_end = rp + pivot_row_length ; 2458 while (rp < rp_end) 2459 { 2460 col = *rp++ ; 2461 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 2462 DEBUG3 (("Col: %d\n", col)) ; 2463 2464 /* clear tags used to construct pivot row pattern */ 2465 col_thickness = -Col [col].shared1.thickness ; 2466 ASSERT (col_thickness > 0) ; 2467 Col [col].shared1.thickness = col_thickness ; 2468 2469 /* === Remove column from degree list =========================== */ 2470 2471 cur_score = Col [col].shared2.score ; 2472 prev_col = Col [col].shared3.prev ; 2473 next_col = Col [col].shared4.degree_next ; 2474 ASSERT (cur_score >= 0) ; 2475 ASSERT (cur_score <= n_col) ; 2476 ASSERT (cur_score >= EMPTY) ; 2477 if (prev_col == EMPTY) 2478 { 2479 head [cur_score] = next_col ; 2480 } 2481 else 2482 { 2483 Col [prev_col].shared4.degree_next = next_col ; 2484 } 2485 if (next_col != EMPTY) 2486 { 2487 Col [next_col].shared3.prev = prev_col ; 2488 } 2489 2490 /* === Scan the column ========================================== */ 2491 2492 cp = &A [Col [col].start] ; 2493 cp_end = cp + Col [col].length ; 2494 while (cp < cp_end) 2495 { 2496 /* get a row */ 2497 row = *cp++ ; 2498 row_mark = Row [row].shared2.mark ; 2499 /* skip if dead */ 2500 if (ROW_IS_MARKED_DEAD (row_mark)) 2501 { 2502 continue ; 2503 } 2504 ASSERT (row != pivot_row) ; 2505 set_difference = row_mark - tag_mark ; 2506 /* check if the row has been seen yet */ 2507 if (set_difference < 0) 2508 { 2509 ASSERT (Row [row].shared1.degree <= max_deg) ; 2510 set_difference = Row [row].shared1.degree ; 2511 } 2512 /* subtract column thickness from this row's set difference */ 2513 set_difference -= col_thickness ; 2514 ASSERT (set_difference >= 0) ; 2515 /* absorb this row if the set difference becomes zero */ 2516 if (set_difference == 0 && aggressive) 2517 { 2518 DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; 2519 KILL_ROW (row) ; 2520 } 2521 else 2522 { 2523 /* save the new mark */ 2524 Row [row].shared2.mark = set_difference + tag_mark ; 2525 } 2526 } 2527 } 2528 2529 #ifndef NDEBUG 2530 debug_deg_lists (n_row, n_col, Row, Col, head, 2531 min_score, n_col2-k-pivot_row_degree, max_deg) ; 2532 #endif /* NDEBUG */ 2533 2534 /* === Add up set differences for each column ======================= */ 2535 2536 DEBUG3 (("** Adding set differences phase. **\n")) ; 2537 2538 /* for each column in pivot row */ 2539 rp = &A [pivot_row_start] ; 2540 rp_end = rp + pivot_row_length ; 2541 while (rp < rp_end) 2542 { 2543 /* get a column */ 2544 col = *rp++ ; 2545 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 2546 hash = 0 ; 2547 cur_score = 0 ; 2548 cp = &A [Col [col].start] ; 2549 /* compact the column */ 2550 new_cp = cp ; 2551 cp_end = cp + Col [col].length ; 2552 2553 DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; 2554 2555 while (cp < cp_end) 2556 { 2557 /* get a row */ 2558 row = *cp++ ; 2559 ASSERT(row >= 0 && row < n_row) ; 2560 row_mark = Row [row].shared2.mark ; 2561 /* skip if dead */ 2562 if (ROW_IS_MARKED_DEAD (row_mark)) 2563 { 2564 DEBUG4 ((" Row %d, dead\n", row)) ; 2565 continue ; 2566 } 2567 DEBUG4 ((" Row %d, set diff %d\n", row, row_mark-tag_mark)); 2568 ASSERT (row_mark >= tag_mark) ; 2569 /* compact the column */ 2570 *new_cp++ = row ; 2571 /* compute hash function */ 2572 hash += row ; 2573 /* add set difference */ 2574 cur_score += row_mark - tag_mark ; 2575 /* integer overflow... */ 2576 cur_score = MIN (cur_score, n_col) ; 2577 } 2578 2579 /* recompute the column's length */ 2580 Col [col].length = (Int) (new_cp - &A [Col [col].start]) ; 2581 2582 /* === Further mass elimination ================================= */ 2583 2584 if (Col [col].length == 0) 2585 { 2586 DEBUG4 (("further mass elimination. Col: %d\n", col)) ; 2587 /* nothing left but the pivot row in this column */ 2588 KILL_PRINCIPAL_COL (col) ; 2589 pivot_row_degree -= Col [col].shared1.thickness ; 2590 ASSERT (pivot_row_degree >= 0) ; 2591 /* order it */ 2592 Col [col].shared2.order = k ; 2593 /* increment order count by column thickness */ 2594 k += Col [col].shared1.thickness ; 2595 } 2596 else 2597 { 2598 /* === Prepare for supercolumn detection ==================== */ 2599 2600 DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; 2601 2602 /* save score so far */ 2603 Col [col].shared2.score = cur_score ; 2604 2605 /* add column to hash table, for supercolumn detection */ 2606 hash %= n_col + 1 ; 2607 2608 DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; 2609 ASSERT (((Int) hash) <= n_col) ; 2610 2611 head_column = head [hash] ; 2612 if (head_column > EMPTY) 2613 { 2614 /* degree list "hash" is non-empty, use prev (shared3) of */ 2615 /* first column in degree list as head of hash bucket */ 2616 first_col = Col [head_column].shared3.headhash ; 2617 Col [head_column].shared3.headhash = col ; 2618 } 2619 else 2620 { 2621 /* degree list "hash" is empty, use head as hash bucket */ 2622 first_col = - (head_column + 2) ; 2623 head [hash] = - (col + 2) ; 2624 } 2625 Col [col].shared4.hash_next = first_col ; 2626 2627 /* save hash function in Col [col].shared3.hash */ 2628 Col [col].shared3.hash = (Int) hash ; 2629 ASSERT (COL_IS_ALIVE (col)) ; 2630 } 2631 } 2632 2633 /* The approximate external column degree is now computed. */ 2634 2635 /* === Supercolumn detection ======================================== */ 2636 2637 DEBUG3 (("** Supercolumn detection phase. **\n")) ; 2638 2639 detect_super_cols ( 2640 2641 #ifndef NDEBUG 2642 n_col, Row, 2643 #endif /* NDEBUG */ 2644 2645 Col, A, head, pivot_row_start, pivot_row_length) ; 2646 2647 /* === Kill the pivotal column ====================================== */ 2648 2649 KILL_PRINCIPAL_COL (pivot_col) ; 2650 2651 /* === Clear mark =================================================== */ 2652 2653 tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ; 2654 2655 #ifndef NDEBUG 2656 DEBUG3 (("check3\n")) ; 2657 debug_mark (n_row, Row, tag_mark, max_mark) ; 2658 #endif /* NDEBUG */ 2659 2660 /* === Finalize the new pivot row, and column scores ================ */ 2661 2662 DEBUG3 (("** Finalize scores phase. **\n")) ; 2663 2664 /* for each column in pivot row */ 2665 rp = &A [pivot_row_start] ; 2666 /* compact the pivot row */ 2667 new_rp = rp ; 2668 rp_end = rp + pivot_row_length ; 2669 while (rp < rp_end) 2670 { 2671 col = *rp++ ; 2672 /* skip dead columns */ 2673 if (COL_IS_DEAD (col)) 2674 { 2675 continue ; 2676 } 2677 *new_rp++ = col ; 2678 /* add new pivot row to column */ 2679 A [Col [col].start + (Col [col].length++)] = pivot_row ; 2680 2681 /* retrieve score so far and add on pivot row's degree. */ 2682 /* (we wait until here for this in case the pivot */ 2683 /* row's degree was reduced due to mass elimination). */ 2684 cur_score = Col [col].shared2.score + pivot_row_degree ; 2685 2686 /* calculate the max possible score as the number of */ 2687 /* external columns minus the 'k' value minus the */ 2688 /* columns thickness */ 2689 max_score = n_col - k - Col [col].shared1.thickness ; 2690 2691 /* make the score the external degree of the union-of-rows */ 2692 cur_score -= Col [col].shared1.thickness ; 2693 2694 /* make sure score is less or equal than the max score */ 2695 cur_score = MIN (cur_score, max_score) ; 2696 ASSERT (cur_score >= 0) ; 2697 2698 /* store updated score */ 2699 Col [col].shared2.score = cur_score ; 2700 2701 /* === Place column back in degree list ========================= */ 2702 2703 ASSERT (min_score >= 0) ; 2704 ASSERT (min_score <= n_col) ; 2705 ASSERT (cur_score >= 0) ; 2706 ASSERT (cur_score <= n_col) ; 2707 ASSERT (head [cur_score] >= EMPTY) ; 2708 next_col = head [cur_score] ; 2709 Col [col].shared4.degree_next = next_col ; 2710 Col [col].shared3.prev = EMPTY ; 2711 if (next_col != EMPTY) 2712 { 2713 Col [next_col].shared3.prev = col ; 2714 } 2715 head [cur_score] = col ; 2716 2717 /* see if this score is less than current min */ 2718 min_score = MIN (min_score, cur_score) ; 2719 2720 } 2721 2722 #ifndef NDEBUG 2723 debug_deg_lists (n_row, n_col, Row, Col, head, 2724 min_score, n_col2-k, max_deg) ; 2725 #endif /* NDEBUG */ 2726 2727 /* === Resurrect the new pivot row ================================== */ 2728 2729 if (pivot_row_degree > 0) 2730 { 2731 /* update pivot row length to reflect any cols that were killed */ 2732 /* during super-col detection and mass elimination */ 2733 Row [pivot_row].start = pivot_row_start ; 2734 Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ; 2735 ASSERT (Row [pivot_row].length > 0) ; 2736 Row [pivot_row].shared1.degree = pivot_row_degree ; 2737 Row [pivot_row].shared2.mark = 0 ; 2738 /* pivot row is no longer dead */ 2739 2740 DEBUG1 (("Resurrect Pivot_row %d deg: %d\n", 2741 pivot_row, pivot_row_degree)) ; 2742 } 2743 } 2744 2745 /* === All principal columns have now been ordered ====================== */ 2746 2747 return (ngarbage) ; 2748 } 2749 2750 2751 /* ========================================================================== */ 2752 /* === order_children ======================================================= */ 2753 /* ========================================================================== */ 2754 2755 /* 2756 The find_ordering routine has ordered all of the principal columns (the 2757 representatives of the supercolumns). The non-principal columns have not 2758 yet been ordered. This routine orders those columns by walking up the 2759 parent tree (a column is a child of the column which absorbed it). The 2760 final permutation vector is then placed in p [0 ... n_col-1], with p [0] 2761 being the first column, and p [n_col-1] being the last. It doesn't look 2762 like it at first glance, but be assured that this routine takes time linear 2763 in the number of columns. Although not immediately obvious, the time 2764 taken by this routine is O (n_col), that is, linear in the number of 2765 columns. Not user-callable. 2766 */ 2767 2768 PRIVATE void order_children 2769 ( 2770 /* === Parameters ======================================================= */ 2771 2772 Int n_col, /* number of columns of A */ 2773 Colamd_Col Col [], /* of size n_col+1 */ 2774 Int p [] /* p [0 ... n_col-1] is the column permutation*/ 2775 ) 2776 { 2777 /* === Local variables ================================================== */ 2778 2779 Int i ; /* loop counter for all columns */ 2780 Int c ; /* column index */ 2781 Int parent ; /* index of column's parent */ 2782 Int order ; /* column's order */ 2783 2784 /* === Order each non-principal column ================================== */ 2785 2786 for (i = 0 ; i < n_col ; i++) 2787 { 2788 /* find an un-ordered non-principal column */ 2789 ASSERT (COL_IS_DEAD (i)) ; 2790 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY) 2791 { 2792 parent = i ; 2793 /* once found, find its principal parent */ 2794 do 2795 { 2796 parent = Col [parent].shared1.parent ; 2797 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; 2798 2799 /* now, order all un-ordered non-principal columns along path */ 2800 /* to this parent. collapse tree at the same time */ 2801 c = i ; 2802 /* get order of parent */ 2803 order = Col [parent].shared2.order ; 2804 2805 do 2806 { 2807 ASSERT (Col [c].shared2.order == EMPTY) ; 2808 2809 /* order this column */ 2810 Col [c].shared2.order = order++ ; 2811 /* collaps tree */ 2812 Col [c].shared1.parent = parent ; 2813 2814 /* get immediate parent of this column */ 2815 c = Col [c].shared1.parent ; 2816 2817 /* continue until we hit an ordered column. There are */ 2818 /* guarranteed not to be anymore unordered columns */ 2819 /* above an ordered column */ 2820 } while (Col [c].shared2.order == EMPTY) ; 2821 2822 /* re-order the super_col parent to largest order for this group */ 2823 Col [parent].shared2.order = order ; 2824 } 2825 } 2826 2827 /* === Generate the permutation ========================================= */ 2828 2829 for (c = 0 ; c < n_col ; c++) 2830 { 2831 p [Col [c].shared2.order] = c ; 2832 } 2833 } 2834 2835 2836 /* ========================================================================== */ 2837 /* === detect_super_cols ==================================================== */ 2838 /* ========================================================================== */ 2839 2840 /* 2841 Detects supercolumns by finding matches between columns in the hash buckets. 2842 Check amongst columns in the set A [row_start ... row_start + row_length-1]. 2843 The columns under consideration are currently *not* in the degree lists, 2844 and have already been placed in the hash buckets. 2845 2846 The hash bucket for columns whose hash function is equal to h is stored 2847 as follows: 2848 2849 if head [h] is >= 0, then head [h] contains a degree list, so: 2850 2851 head [h] is the first column in degree bucket h. 2852 Col [head [h]].headhash gives the first column in hash bucket h. 2853 2854 otherwise, the degree list is empty, and: 2855 2856 -(head [h] + 2) is the first column in hash bucket h. 2857 2858 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous 2859 column" pointer. Col [c].shared3.hash is used instead as the hash number 2860 for that column. The value of Col [c].shared4.hash_next is the next column 2861 in the same hash bucket. 2862 2863 Assuming no, or "few" hash collisions, the time taken by this routine is 2864 linear in the sum of the sizes (lengths) of each column whose score has 2865 just been computed in the approximate degree computation. 2866 Not user-callable. 2867 */ 2868 2869 PRIVATE void detect_super_cols 2870 ( 2871 /* === Parameters ======================================================= */ 2872 2873 #ifndef NDEBUG 2874 /* these two parameters are only needed when debugging is enabled: */ 2875 Int n_col, /* number of columns of A */ 2876 Colamd_Row Row [], /* of size n_row+1 */ 2877 #endif /* NDEBUG */ 2878 2879 Colamd_Col Col [], /* of size n_col+1 */ 2880 Int A [], /* row indices of A */ 2881 Int head [], /* head of degree lists and hash buckets */ 2882 Int row_start, /* pointer to set of columns to check */ 2883 Int row_length /* number of columns to check */ 2884 ) 2885 { 2886 /* === Local variables ================================================== */ 2887 2888 Int hash ; /* hash value for a column */ 2889 Int *rp ; /* pointer to a row */ 2890 Int c ; /* a column index */ 2891 Int super_c ; /* column index of the column to absorb into */ 2892 Int *cp1 ; /* column pointer for column super_c */ 2893 Int *cp2 ; /* column pointer for column c */ 2894 Int length ; /* length of column super_c */ 2895 Int prev_c ; /* column preceding c in hash bucket */ 2896 Int i ; /* loop counter */ 2897 Int *rp_end ; /* pointer to the end of the row */ 2898 Int col ; /* a column index in the row to check */ 2899 Int head_column ; /* first column in hash bucket or degree list */ 2900 Int first_col ; /* first column in hash bucket */ 2901 2902 /* === Consider each column in the row ================================== */ 2903 2904 rp = &A [row_start] ; 2905 rp_end = rp + row_length ; 2906 while (rp < rp_end) 2907 { 2908 col = *rp++ ; 2909 if (COL_IS_DEAD (col)) 2910 { 2911 continue ; 2912 } 2913 2914 /* get hash number for this column */ 2915 hash = Col [col].shared3.hash ; 2916 ASSERT (hash <= n_col) ; 2917 2918 /* === Get the first column in this hash bucket ===================== */ 2919 2920 head_column = head [hash] ; 2921 if (head_column > EMPTY) 2922 { 2923 first_col = Col [head_column].shared3.headhash ; 2924 } 2925 else 2926 { 2927 first_col = - (head_column + 2) ; 2928 } 2929 2930 /* === Consider each column in the hash bucket ====================== */ 2931 2932 for (super_c = first_col ; super_c != EMPTY ; 2933 super_c = Col [super_c].shared4.hash_next) 2934 { 2935 ASSERT (COL_IS_ALIVE (super_c)) ; 2936 ASSERT (Col [super_c].shared3.hash == hash) ; 2937 length = Col [super_c].length ; 2938 2939 /* prev_c is the column preceding column c in the hash bucket */ 2940 prev_c = super_c ; 2941 2942 /* === Compare super_c with all columns after it ================ */ 2943 2944 for (c = Col [super_c].shared4.hash_next ; 2945 c != EMPTY ; c = Col [c].shared4.hash_next) 2946 { 2947 ASSERT (c != super_c) ; 2948 ASSERT (COL_IS_ALIVE (c)) ; 2949 ASSERT (Col [c].shared3.hash == hash) ; 2950 2951 /* not identical if lengths or scores are different */ 2952 if (Col [c].length != length || 2953 Col [c].shared2.score != Col [super_c].shared2.score) 2954 { 2955 prev_c = c ; 2956 continue ; 2957 } 2958 2959 /* compare the two columns */ 2960 cp1 = &A [Col [super_c].start] ; 2961 cp2 = &A [Col [c].start] ; 2962 2963 for (i = 0 ; i < length ; i++) 2964 { 2965 /* the columns are "clean" (no dead rows) */ 2966 ASSERT (ROW_IS_ALIVE (*cp1)) ; 2967 ASSERT (ROW_IS_ALIVE (*cp2)) ; 2968 /* row indices will same order for both supercols, */ 2969 /* no gather scatter nessasary */ 2970 if (*cp1++ != *cp2++) 2971 { 2972 break ; 2973 } 2974 } 2975 2976 /* the two columns are different if the for-loop "broke" */ 2977 if (i != length) 2978 { 2979 prev_c = c ; 2980 continue ; 2981 } 2982 2983 /* === Got it! two columns are identical =================== */ 2984 2985 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; 2986 2987 Col [super_c].shared1.thickness += Col [c].shared1.thickness ; 2988 Col [c].shared1.parent = super_c ; 2989 KILL_NON_PRINCIPAL_COL (c) ; 2990 /* order c later, in order_children() */ 2991 Col [c].shared2.order = EMPTY ; 2992 /* remove c from hash bucket */ 2993 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; 2994 } 2995 } 2996 2997 /* === Empty this hash bucket ======================================= */ 2998 2999 if (head_column > EMPTY) 3000 { 3001 /* corresponding degree list "hash" is not empty */ 3002 Col [head_column].shared3.headhash = EMPTY ; 3003 } 3004 else 3005 { 3006 /* corresponding degree list "hash" is empty */ 3007 head [hash] = EMPTY ; 3008 } 3009 } 3010 } 3011 3012 3013 /* ========================================================================== */ 3014 /* === garbage_collection =================================================== */ 3015 /* ========================================================================== */ 3016 3017 /* 3018 Defragments and compacts columns and rows in the workspace A. Used when 3019 all avaliable memory has been used while performing row merging. Returns 3020 the index of the first free position in A, after garbage collection. The 3021 time taken by this routine is linear is the size of the array A, which is 3022 itself linear in the number of nonzeros in the input matrix. 3023 Not user-callable. 3024 */ 3025 3026 PRIVATE Int garbage_collection /* returns the new value of pfree */ 3027 ( 3028 /* === Parameters ======================================================= */ 3029 3030 Int n_row, /* number of rows */ 3031 Int n_col, /* number of columns */ 3032 Colamd_Row Row [], /* row info */ 3033 Colamd_Col Col [], /* column info */ 3034 Int A [], /* A [0 ... Alen-1] holds the matrix */ 3035 Int *pfree /* &A [0] ... pfree is in use */ 3036 ) 3037 { 3038 /* === Local variables ================================================== */ 3039 3040 Int *psrc ; /* source pointer */ 3041 Int *pdest ; /* destination pointer */ 3042 Int j ; /* counter */ 3043 Int r ; /* a row index */ 3044 Int c ; /* a column index */ 3045 Int length ; /* length of a row or column */ 3046 3047 #ifndef NDEBUG 3048 Int debug_rows ; 3049 DEBUG2 (("Defrag..\n")) ; 3050 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ; 3051 debug_rows = 0 ; 3052 #endif /* NDEBUG */ 3053 3054 /* === Defragment the columns =========================================== */ 3055 3056 pdest = &A[0] ; 3057 for (c = 0 ; c < n_col ; c++) 3058 { 3059 if (COL_IS_ALIVE (c)) 3060 { 3061 psrc = &A [Col [c].start] ; 3062 3063 /* move and compact the column */ 3064 ASSERT (pdest <= psrc) ; 3065 Col [c].start = (Int) (pdest - &A [0]) ; 3066 length = Col [c].length ; 3067 for (j = 0 ; j < length ; j++) 3068 { 3069 r = *psrc++ ; 3070 if (ROW_IS_ALIVE (r)) 3071 { 3072 *pdest++ = r ; 3073 } 3074 } 3075 Col [c].length = (Int) (pdest - &A [Col [c].start]) ; 3076 } 3077 } 3078 3079 /* === Prepare to defragment the rows =================================== */ 3080 3081 for (r = 0 ; r < n_row ; r++) 3082 { 3083 if (ROW_IS_DEAD (r) || (Row [r].length == 0)) 3084 { 3085 /* This row is already dead, or is of zero length. Cannot compact 3086 * a row of zero length, so kill it. NOTE: in the current version, 3087 * there are no zero-length live rows. Kill the row (for the first 3088 * time, or again) just to be safe. */ 3089 KILL_ROW (r) ; 3090 } 3091 else 3092 { 3093 /* save first column index in Row [r].shared2.first_column */ 3094 psrc = &A [Row [r].start] ; 3095 Row [r].shared2.first_column = *psrc ; 3096 ASSERT (ROW_IS_ALIVE (r)) ; 3097 /* flag the start of the row with the one's complement of row */ 3098 *psrc = ONES_COMPLEMENT (r) ; 3099 #ifndef NDEBUG 3100 debug_rows++ ; 3101 #endif /* NDEBUG */ 3102 } 3103 } 3104 3105 /* === Defragment the rows ============================================== */ 3106 3107 psrc = pdest ; 3108 while (psrc < pfree) 3109 { 3110 /* find a negative number ... the start of a row */ 3111 if (*psrc++ < 0) 3112 { 3113 psrc-- ; 3114 /* get the row index */ 3115 r = ONES_COMPLEMENT (*psrc) ; 3116 ASSERT (r >= 0 && r < n_row) ; 3117 /* restore first column index */ 3118 *psrc = Row [r].shared2.first_column ; 3119 ASSERT (ROW_IS_ALIVE (r)) ; 3120 ASSERT (Row [r].length > 0) ; 3121 /* move and compact the row */ 3122 ASSERT (pdest <= psrc) ; 3123 Row [r].start = (Int) (pdest - &A [0]) ; 3124 length = Row [r].length ; 3125 for (j = 0 ; j < length ; j++) 3126 { 3127 c = *psrc++ ; 3128 if (COL_IS_ALIVE (c)) 3129 { 3130 *pdest++ = c ; 3131 } 3132 } 3133 Row [r].length = (Int) (pdest - &A [Row [r].start]) ; 3134 ASSERT (Row [r].length > 0) ; 3135 #ifndef NDEBUG 3136 debug_rows-- ; 3137 #endif /* NDEBUG */ 3138 } 3139 } 3140 /* ensure we found all the rows */ 3141 ASSERT (debug_rows == 0) ; 3142 3143 /* === Return the new value of pfree ==================================== */ 3144 3145 return ((Int) (pdest - &A [0])) ; 3146 } 3147 3148 3149 /* ========================================================================== */ 3150 /* === clear_mark =========================================================== */ 3151 /* ========================================================================== */ 3152 3153 /* 3154 Clears the Row [].shared2.mark array, and returns the new tag_mark. 3155 Return value is the new tag_mark. Not user-callable. 3156 */ 3157 3158 PRIVATE Int clear_mark /* return the new value for tag_mark */ 3159 ( 3160 /* === Parameters ======================================================= */ 3161 3162 Int tag_mark, /* new value of tag_mark */ 3163 Int max_mark, /* max allowed value of tag_mark */ 3164 3165 Int n_row, /* number of rows in A */ 3166 Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ 3167 ) 3168 { 3169 /* === Local variables ================================================== */ 3170 3171 Int r ; 3172 3173 if (tag_mark <= 0 || tag_mark >= max_mark) 3174 { 3175 for (r = 0 ; r < n_row ; r++) 3176 { 3177 if (ROW_IS_ALIVE (r)) 3178 { 3179 Row [r].shared2.mark = 0 ; 3180 } 3181 } 3182 tag_mark = 1 ; 3183 } 3184 3185 return (tag_mark) ; 3186 } 3187 3188 3189 /* ========================================================================== */ 3190 /* === print_report ========================================================= */ 3191 /* ========================================================================== */ 3192 3193 PRIVATE void print_report 3194 ( 3195 char *method, 3196 Int stats [COLAMD_STATS] 3197 ) 3198 { 3199 3200 Int i1, i2, i3 ; 3201 3202 PRINTF (("\n%s version %d.%d, %s: ", method, 3203 COLAMD_MAIN_VERSION, COLAMD_SUB_VERSION, COLAMD_DATE)) ; 3204 3205 if (!stats) 3206 { 3207 PRINTF (("No statistics available.\n")) ; 3208 return ; 3209 } 3210 3211 i1 = stats [COLAMD_INFO1] ; 3212 i2 = stats [COLAMD_INFO2] ; 3213 i3 = stats [COLAMD_INFO3] ; 3214 3215 if (stats [COLAMD_STATUS] >= 0) 3216 { 3217 PRINTF (("OK. ")) ; 3218 } 3219 else 3220 { 3221 PRINTF (("ERROR. ")) ; 3222 } 3223 3224 switch (stats [COLAMD_STATUS]) 3225 { 3226 3227 case COLAMD_OK_BUT_JUMBLED: 3228 3229 PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ; 3230 3231 PRINTF(("%s: number of duplicate or out-of-order row indices: %d\n", 3232 method, i3)) ; 3233 3234 PRINTF(("%s: last seen duplicate or out-of-order row index: %d\n", 3235 method, INDEX (i2))) ; 3236 3237 PRINTF(("%s: last seen in column: %d", 3238 method, INDEX (i1))) ; 3239 3240 /* no break - fall through to next case instead */ 3241 3242 case COLAMD_OK: 3243 3244 PRINTF(("\n")) ; 3245 3246 PRINTF(("%s: number of dense or empty rows ignored: %d\n", 3247 method, stats [COLAMD_DENSE_ROW])) ; 3248 3249 PRINTF(("%s: number of dense or empty columns ignored: %d\n", 3250 method, stats [COLAMD_DENSE_COL])) ; 3251 3252 PRINTF(("%s: number of garbage collections performed: %d\n", 3253 method, stats [COLAMD_DEFRAG_COUNT])) ; 3254 break ; 3255 3256 case COLAMD_ERROR_A_not_present: 3257 3258 PRINTF(("Array A (row indices of matrix) not present.\n")) ; 3259 break ; 3260 3261 case COLAMD_ERROR_p_not_present: 3262 3263 PRINTF(("Array p (column pointers for matrix) not present.\n")) ; 3264 break ; 3265 3266 case COLAMD_ERROR_nrow_negative: 3267 3268 PRINTF(("Invalid number of rows (%d).\n", i1)) ; 3269 break ; 3270 3271 case COLAMD_ERROR_ncol_negative: 3272 3273 PRINTF(("Invalid number of columns (%d).\n", i1)) ; 3274 break ; 3275 3276 case COLAMD_ERROR_nnz_negative: 3277 3278 PRINTF(("Invalid number of nonzero entries (%d).\n", i1)) ; 3279 break ; 3280 3281 case COLAMD_ERROR_p0_nonzero: 3282 3283 PRINTF(("Invalid column pointer, p [0] = %d, must be zero.\n", i1)); 3284 break ; 3285 3286 case COLAMD_ERROR_A_too_small: 3287 3288 PRINTF(("Array A too small.\n")) ; 3289 PRINTF((" Need Alen >= %d, but given only Alen = %d.\n", 3290 i1, i2)) ; 3291 break ; 3292 3293 case COLAMD_ERROR_col_length_negative: 3294 3295 PRINTF 3296 (("Column %d has a negative number of nonzero entries (%d).\n", 3297 INDEX (i1), i2)) ; 3298 break ; 3299 3300 case COLAMD_ERROR_row_index_out_of_bounds: 3301 3302 PRINTF 3303 (("Row index (row %d) out of bounds (%d to %d) in column %d.\n", 3304 INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1))) ; 3305 break ; 3306 3307 case COLAMD_ERROR_out_of_memory: 3308 3309 PRINTF(("Out of memory.\n")) ; 3310 break ; 3311 3312 /* v2.4: internal-error case deleted */ 3313 } 3314 } 3315 3316 3317 3318 3319 /* ========================================================================== */ 3320 /* === colamd debugging routines ============================================ */ 3321 /* ========================================================================== */ 3322 3323 /* When debugging is disabled, the remainder of this file is ignored. */ 3324 3325 #ifndef NDEBUG 3326 3327 3328 /* ========================================================================== */ 3329 /* === debug_structures ===================================================== */ 3330 /* ========================================================================== */ 3331 3332 /* 3333 At this point, all empty rows and columns are dead. All live columns 3334 are "clean" (containing no dead rows) and simplicial (no supercolumns 3335 yet). Rows may contain dead columns, but all live rows contain at 3336 least one live column. 3337 */ 3338 3339 PRIVATE void debug_structures 3340 ( 3341 /* === Parameters ======================================================= */ 3342 3343 Int n_row, 3344 Int n_col, 3345 Colamd_Row Row [], 3346 Colamd_Col Col [], 3347 Int A [], 3348 Int n_col2 3349 ) 3350 { 3351 /* === Local variables ================================================== */ 3352 3353 Int i ; 3354 Int c ; 3355 Int *cp ; 3356 Int *cp_end ; 3357 Int len ; 3358 Int score ; 3359 Int r ; 3360 Int *rp ; 3361 Int *rp_end ; 3362 Int deg ; 3363 3364 /* === Check A, Row, and Col ============================================ */ 3365 3366 for (c = 0 ; c < n_col ; c++) 3367 { 3368 if (COL_IS_ALIVE (c)) 3369 { 3370 len = Col [c].length ; 3371 score = Col [c].shared2.score ; 3372 DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; 3373 ASSERT (len > 0) ; 3374 ASSERT (score >= 0) ; 3375 ASSERT (Col [c].shared1.thickness == 1) ; 3376 cp = &A [Col [c].start] ; 3377 cp_end = cp + len ; 3378 while (cp < cp_end) 3379 { 3380 r = *cp++ ; 3381 ASSERT (ROW_IS_ALIVE (r)) ; 3382 } 3383 } 3384 else 3385 { 3386 i = Col [c].shared2.order ; 3387 ASSERT (i >= n_col2 && i < n_col) ; 3388 } 3389 } 3390 3391 for (r = 0 ; r < n_row ; r++) 3392 { 3393 if (ROW_IS_ALIVE (r)) 3394 { 3395 i = 0 ; 3396 len = Row [r].length ; 3397 deg = Row [r].shared1.degree ; 3398 ASSERT (len > 0) ; 3399 ASSERT (deg > 0) ; 3400 rp = &A [Row [r].start] ; 3401 rp_end = rp + len ; 3402 while (rp < rp_end) 3403 { 3404 c = *rp++ ; 3405 if (COL_IS_ALIVE (c)) 3406 { 3407 i++ ; 3408 } 3409 } 3410 ASSERT (i > 0) ; 3411 } 3412 } 3413 } 3414 3415 3416 /* ========================================================================== */ 3417 /* === debug_deg_lists ====================================================== */ 3418 /* ========================================================================== */ 3419 3420 /* 3421 Prints the contents of the degree lists. Counts the number of columns 3422 in the degree list and compares it to the total it should have. Also 3423 checks the row degrees. 3424 */ 3425 3426 PRIVATE void debug_deg_lists 3427 ( 3428 /* === Parameters ======================================================= */ 3429 3430 Int n_row, 3431 Int n_col, 3432 Colamd_Row Row [], 3433 Colamd_Col Col [], 3434 Int head [], 3435 Int min_score, 3436 Int should, 3437 Int max_deg 3438 ) 3439 { 3440 /* === Local variables ================================================== */ 3441 3442 Int deg ; 3443 Int col ; 3444 Int have ; 3445 Int row ; 3446 3447 /* === Check the degree lists =========================================== */ 3448 3449 if (n_col > 10000 && colamd_debug <= 0) 3450 { 3451 return ; 3452 } 3453 have = 0 ; 3454 DEBUG4 (("Degree lists: %d\n", min_score)) ; 3455 for (deg = 0 ; deg <= n_col ; deg++) 3456 { 3457 col = head [deg] ; 3458 if (col == EMPTY) 3459 { 3460 continue ; 3461 } 3462 DEBUG4 (("%d:", deg)) ; 3463 while (col != EMPTY) 3464 { 3465 DEBUG4 ((" %d", col)) ; 3466 have += Col [col].shared1.thickness ; 3467 ASSERT (COL_IS_ALIVE (col)) ; 3468 col = Col [col].shared4.degree_next ; 3469 } 3470 DEBUG4 (("\n")) ; 3471 } 3472 DEBUG4 (("should %d have %d\n", should, have)) ; 3473 ASSERT (should == have) ; 3474 3475 /* === Check the row degrees ============================================ */ 3476 3477 if (n_row > 10000 && colamd_debug <= 0) 3478 { 3479 return ; 3480 } 3481 for (row = 0 ; row < n_row ; row++) 3482 { 3483 if (ROW_IS_ALIVE (row)) 3484 { 3485 ASSERT (Row [row].shared1.degree <= max_deg) ; 3486 } 3487 } 3488 } 3489 3490 3491 /* ========================================================================== */ 3492 /* === debug_mark =========================================================== */ 3493 /* ========================================================================== */ 3494 3495 /* 3496 Ensures that the tag_mark is less that the maximum and also ensures that 3497 each entry in the mark array is less than the tag mark. 3498 */ 3499 3500 PRIVATE void debug_mark 3501 ( 3502 /* === Parameters ======================================================= */ 3503 3504 Int n_row, 3505 Colamd_Row Row [], 3506 Int tag_mark, 3507 Int max_mark 3508 ) 3509 { 3510 /* === Local variables ================================================== */ 3511 3512 Int r ; 3513 3514 /* === Check the Row marks ============================================== */ 3515 3516 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; 3517 if (n_row > 10000 && colamd_debug <= 0) 3518 { 3519 return ; 3520 } 3521 for (r = 0 ; r < n_row ; r++) 3522 { 3523 ASSERT (Row [r].shared2.mark < tag_mark) ; 3524 } 3525 } 3526 3527 3528 /* ========================================================================== */ 3529 /* === debug_matrix ========================================================= */ 3530 /* ========================================================================== */ 3531 3532 /* 3533 Prints out the contents of the columns and the rows. 3534 */ 3535 3536 PRIVATE void debug_matrix 3537 ( 3538 /* === Parameters ======================================================= */ 3539 3540 Int n_row, 3541 Int n_col, 3542 Colamd_Row Row [], 3543 Colamd_Col Col [], 3544 Int A [] 3545 ) 3546 { 3547 /* === Local variables ================================================== */ 3548 3549 Int r ; 3550 Int c ; 3551 Int *rp ; 3552 Int *rp_end ; 3553 Int *cp ; 3554 Int *cp_end ; 3555 3556 /* === Dump the rows and columns of the matrix ========================== */ 3557 3558 if (colamd_debug < 3) 3559 { 3560 return ; 3561 } 3562 DEBUG3 (("DUMP MATRIX:\n")) ; 3563 for (r = 0 ; r < n_row ; r++) 3564 { 3565 DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ; 3566 if (ROW_IS_DEAD (r)) 3567 { 3568 continue ; 3569 } 3570 DEBUG3 (("start %d length %d degree %d\n", 3571 Row [r].start, Row [r].length, Row [r].shared1.degree)) ; 3572 rp = &A [Row [r].start] ; 3573 rp_end = rp + Row [r].length ; 3574 while (rp < rp_end) 3575 { 3576 c = *rp++ ; 3577 DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ; 3578 } 3579 } 3580 3581 for (c = 0 ; c < n_col ; c++) 3582 { 3583 DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ; 3584 if (COL_IS_DEAD (c)) 3585 { 3586 continue ; 3587 } 3588 DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", 3589 Col [c].start, Col [c].length, 3590 Col [c].shared1.thickness, Col [c].shared2.score)) ; 3591 cp = &A [Col [c].start] ; 3592 cp_end = cp + Col [c].length ; 3593 while (cp < cp_end) 3594 { 3595 r = *cp++ ; 3596 DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ; 3597 } 3598 } 3599 } 3600 3601 PRIVATE void colamd_get_debug 3602 ( 3603 char *method 3604 ) 3605 { 3606 FILE *f ; 3607 colamd_debug = 0 ; /* no debug printing */ 3608 f = fopen ("debug", "r") ; 3609 if (f == (FILE *) NULL) 3610 { 3611 colamd_debug = 0 ; 3612 } 3613 else 3614 { 3615 fscanf (f, "%d", &colamd_debug) ; 3616 fclose (f) ; 3617 } 3618 DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", 3619 method, colamd_debug)) ; 3620 } 3621 3622 #endif /* NDEBUG */ 3623