1// // This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> 5// 6// This Source Code Form is subject to the terms of the Mozilla 7// Public License v. 2.0. If a copy of the MPL was not distributed 8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10// This file is modified from the colamd/symamd library. The copyright is below 11 12// The authors of the code itself are Stefan I. Larimore and Timothy A. 13// Davis (davis@cise.ufl.edu), University of Florida. The algorithm was 14// developed in collaboration with John Gilbert, Xerox PARC, and Esmond 15// Ng, Oak Ridge National Laboratory. 16// 17// Date: 18// 19// September 8, 2003. Version 2.3. 20// 21// Acknowledgements: 22// 23// This work was supported by the National Science Foundation, under 24// grants DMS-9504974 and DMS-9803599. 25// 26// Notice: 27// 28// Copyright (c) 1998-2003 by the University of Florida. 29// All Rights Reserved. 30// 31// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 32// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 33// 34// Permission is hereby granted to use, copy, modify, and/or distribute 35// this program, provided that the Copyright, this License, and the 36// Availability of the original version is retained on all copies and made 37// accessible to the end-user of any code or package that includes COLAMD 38// or any modified version of COLAMD. 39// 40// Availability: 41// 42// The colamd/symamd library is available at 43// 44// http://www.suitesparse.com 45 46 47#ifndef EIGEN_COLAMD_H 48#define EIGEN_COLAMD_H 49 50namespace internal { 51/* Ensure that debugging is turned off: */ 52#ifndef COLAMD_NDEBUG 53#define COLAMD_NDEBUG 54#endif /* NDEBUG */ 55/* ========================================================================== */ 56/* === Knob and statistics definitions ====================================== */ 57/* ========================================================================== */ 58 59/* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ 60#define COLAMD_KNOBS 20 61 62/* number of output statistics. Only stats [0..6] are currently used. */ 63#define COLAMD_STATS 20 64 65/* knobs [0] and stats [0]: dense row knob and output statistic. */ 66#define COLAMD_DENSE_ROW 0 67 68/* knobs [1] and stats [1]: dense column knob and output statistic. */ 69#define COLAMD_DENSE_COL 1 70 71/* stats [2]: memory defragmentation count output statistic */ 72#define COLAMD_DEFRAG_COUNT 2 73 74/* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ 75#define COLAMD_STATUS 3 76 77/* stats [4..6]: error info, or info on jumbled columns */ 78#define COLAMD_INFO1 4 79#define COLAMD_INFO2 5 80#define COLAMD_INFO3 6 81 82/* error codes returned in stats [3]: */ 83#define COLAMD_OK (0) 84#define COLAMD_OK_BUT_JUMBLED (1) 85#define COLAMD_ERROR_A_not_present (-1) 86#define COLAMD_ERROR_p_not_present (-2) 87#define COLAMD_ERROR_nrow_negative (-3) 88#define COLAMD_ERROR_ncol_negative (-4) 89#define COLAMD_ERROR_nnz_negative (-5) 90#define COLAMD_ERROR_p0_nonzero (-6) 91#define COLAMD_ERROR_A_too_small (-7) 92#define COLAMD_ERROR_col_length_negative (-8) 93#define COLAMD_ERROR_row_index_out_of_bounds (-9) 94#define COLAMD_ERROR_out_of_memory (-10) 95#define COLAMD_ERROR_internal_error (-999) 96 97/* ========================================================================== */ 98/* === Definitions ========================================================== */ 99/* ========================================================================== */ 100 101#define ONES_COMPLEMENT(r) (-(r)-1) 102 103/* -------------------------------------------------------------------------- */ 104 105#define COLAMD_EMPTY (-1) 106 107/* Row and column status */ 108#define ALIVE (0) 109#define DEAD (-1) 110 111/* Column status */ 112#define DEAD_PRINCIPAL (-1) 113#define DEAD_NON_PRINCIPAL (-2) 114 115/* Macros for row and column status update and checking. */ 116#define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) 117#define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) 118#define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) 119#define COL_IS_DEAD(c) (Col [c].start < ALIVE) 120#define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) 121#define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) 122#define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } 123#define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } 124#define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } 125 126/* ========================================================================== */ 127/* === Colamd reporting mechanism =========================================== */ 128/* ========================================================================== */ 129 130// == Row and Column structures == 131template <typename IndexType> 132struct colamd_col 133{ 134 IndexType start ; /* index for A of first row in this column, or DEAD */ 135 /* if column is dead */ 136 IndexType length ; /* number of rows in this column */ 137 union 138 { 139 IndexType thickness ; /* number of original columns represented by this */ 140 /* col, if the column is alive */ 141 IndexType parent ; /* parent in parent tree super-column structure, if */ 142 /* the column is dead */ 143 } shared1 ; 144 union 145 { 146 IndexType score ; /* the score used to maintain heap, if col is alive */ 147 IndexType order ; /* pivot ordering of this column, if col is dead */ 148 } shared2 ; 149 union 150 { 151 IndexType headhash ; /* head of a hash bucket, if col is at the head of */ 152 /* a degree list */ 153 IndexType hash ; /* hash value, if col is not in a degree list */ 154 IndexType prev ; /* previous column in degree list, if col is in a */ 155 /* degree list (but not at the head of a degree list) */ 156 } shared3 ; 157 union 158 { 159 IndexType degree_next ; /* next column, if col is in a degree list */ 160 IndexType hash_next ; /* next column, if col is in a hash list */ 161 } shared4 ; 162 163}; 164 165template <typename IndexType> 166struct Colamd_Row 167{ 168 IndexType start ; /* index for A of first col in this row */ 169 IndexType length ; /* number of principal columns in this row */ 170 union 171 { 172 IndexType degree ; /* number of principal & non-principal columns in row */ 173 IndexType p ; /* used as a row pointer in init_rows_cols () */ 174 } shared1 ; 175 union 176 { 177 IndexType mark ; /* for computing set differences and marking dead rows*/ 178 IndexType first_column ;/* first column in row (used in garbage collection) */ 179 } shared2 ; 180 181}; 182 183/* ========================================================================== */ 184/* === Colamd recommended memory size ======================================= */ 185/* ========================================================================== */ 186 187/* 188 The recommended length Alen of the array A passed to colamd is given by 189 the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any 190 argument is negative. 2*nnz space is required for the row and column 191 indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is 192 required for the Col and Row arrays, respectively, which are internal to 193 colamd. An additional n_col space is the minimal amount of "elbow room", 194 and nnz/5 more space is recommended for run time efficiency. 195 196 This macro is not needed when using symamd. 197 198 Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid 199 gcc -pedantic warning messages. 200*/ 201template <typename IndexType> 202inline IndexType colamd_c(IndexType n_col) 203{ return IndexType( ((n_col) + 1) * sizeof (colamd_col<IndexType>) / sizeof (IndexType) ) ; } 204 205template <typename IndexType> 206inline IndexType colamd_r(IndexType n_row) 207{ return IndexType(((n_row) + 1) * sizeof (Colamd_Row<IndexType>) / sizeof (IndexType)); } 208 209// Prototypes of non-user callable routines 210template <typename IndexType> 211static IndexType init_rows_cols (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> col [], IndexType A [], IndexType p [], IndexType stats[COLAMD_STATS] ); 212 213template <typename IndexType> 214static void init_scoring (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], double knobs[COLAMD_KNOBS], IndexType *p_n_row2, IndexType *p_n_col2, IndexType *p_max_deg); 215 216template <typename IndexType> 217static IndexType find_ordering (IndexType n_row, IndexType n_col, IndexType Alen, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType n_col2, IndexType max_deg, IndexType pfree); 218 219template <typename IndexType> 220static void order_children (IndexType n_col, colamd_col<IndexType> Col [], IndexType p []); 221 222template <typename IndexType> 223static void detect_super_cols (colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType row_start, IndexType row_length ) ; 224 225template <typename IndexType> 226static IndexType garbage_collection (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType *pfree) ; 227 228template <typename IndexType> 229static inline IndexType clear_mark (IndexType n_row, Colamd_Row<IndexType> Row [] ) ; 230 231/* === No debugging ========================================================= */ 232 233#define COLAMD_DEBUG0(params) ; 234#define COLAMD_DEBUG1(params) ; 235#define COLAMD_DEBUG2(params) ; 236#define COLAMD_DEBUG3(params) ; 237#define COLAMD_DEBUG4(params) ; 238 239#define COLAMD_ASSERT(expression) ((void) 0) 240 241 242/** 243 * \brief Returns the recommended value of Alen 244 * 245 * Returns recommended value of Alen for use by colamd. 246 * Returns -1 if any input argument is negative. 247 * The use of this routine or macro is optional. 248 * Note that the macro uses its arguments more than once, 249 * so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED. 250 * 251 * \param nnz nonzeros in A 252 * \param n_row number of rows in A 253 * \param n_col number of columns in A 254 * \return recommended value of Alen for use by colamd 255 */ 256template <typename IndexType> 257inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col) 258{ 259 if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) 260 return (-1); 261 else 262 return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5)); 263} 264 265/** 266 * \brief set default parameters The use of this routine is optional. 267 * 268 * Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) 269 * entries are removed prior to ordering. Columns with more than 270 * (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to 271 * ordering, and placed last in the output column ordering. 272 * 273 * COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1, 274 * respectively, in colamd.h. Default values of these two knobs 275 * are both 0.5. Currently, only knobs [0] and knobs [1] are 276 * used, but future versions may use more knobs. If so, they will 277 * be properly set to their defaults by the future version of 278 * colamd_set_defaults, so that the code that calls colamd will 279 * not need to change, assuming that you either use 280 * colamd_set_defaults, or pass a (double *) NULL pointer as the 281 * knobs array to colamd or symamd. 282 * 283 * \param knobs parameter settings for colamd 284 */ 285 286static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS]) 287{ 288 /* === Local variables ================================================== */ 289 290 int i ; 291 292 if (!knobs) 293 { 294 return ; /* no knobs to initialize */ 295 } 296 for (i = 0 ; i < COLAMD_KNOBS ; i++) 297 { 298 knobs [i] = 0 ; 299 } 300 knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ 301 knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ 302} 303 304/** 305 * \brief Computes a column ordering using the column approximate minimum degree ordering 306 * 307 * Computes a column ordering (Q) of A such that P(AQ)=LU or 308 * (AQ)'AQ=LL' have less fill-in and require fewer floating point 309 * operations than factorizing the unpermuted matrix A or A'A, 310 * respectively. 311 * 312 * 313 * \param n_row number of rows in A 314 * \param n_col number of columns in A 315 * \param Alen, size of the array A 316 * \param A row indices of the matrix, of size ALen 317 * \param p column pointers of A, of size n_col+1 318 * \param knobs parameter settings for colamd 319 * \param stats colamd output statistics and error codes 320 */ 321template <typename IndexType> 322static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS]) 323{ 324 /* === Local variables ================================================== */ 325 326 IndexType i ; /* loop index */ 327 IndexType nnz ; /* nonzeros in A */ 328 IndexType Row_size ; /* size of Row [], in integers */ 329 IndexType Col_size ; /* size of Col [], in integers */ 330 IndexType need ; /* minimum required length of A */ 331 Colamd_Row<IndexType> *Row ; /* pointer into A of Row [0..n_row] array */ 332 colamd_col<IndexType> *Col ; /* pointer into A of Col [0..n_col] array */ 333 IndexType n_col2 ; /* number of non-dense, non-empty columns */ 334 IndexType n_row2 ; /* number of non-dense, non-empty rows */ 335 IndexType ngarbage ; /* number of garbage collections performed */ 336 IndexType max_deg ; /* maximum row degree */ 337 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ 338 339 340 /* === Check the input arguments ======================================== */ 341 342 if (!stats) 343 { 344 COLAMD_DEBUG0 (("colamd: stats not present\n")) ; 345 return (false) ; 346 } 347 for (i = 0 ; i < COLAMD_STATS ; i++) 348 { 349 stats [i] = 0 ; 350 } 351 stats [COLAMD_STATUS] = COLAMD_OK ; 352 stats [COLAMD_INFO1] = -1 ; 353 stats [COLAMD_INFO2] = -1 ; 354 355 if (!A) /* A is not present */ 356 { 357 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; 358 COLAMD_DEBUG0 (("colamd: A not present\n")) ; 359 return (false) ; 360 } 361 362 if (!p) /* p is not present */ 363 { 364 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; 365 COLAMD_DEBUG0 (("colamd: p not present\n")) ; 366 return (false) ; 367 } 368 369 if (n_row < 0) /* n_row must be >= 0 */ 370 { 371 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; 372 stats [COLAMD_INFO1] = n_row ; 373 COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; 374 return (false) ; 375 } 376 377 if (n_col < 0) /* n_col must be >= 0 */ 378 { 379 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; 380 stats [COLAMD_INFO1] = n_col ; 381 COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; 382 return (false) ; 383 } 384 385 nnz = p [n_col] ; 386 if (nnz < 0) /* nnz must be >= 0 */ 387 { 388 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; 389 stats [COLAMD_INFO1] = nnz ; 390 COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; 391 return (false) ; 392 } 393 394 if (p [0] != 0) 395 { 396 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; 397 stats [COLAMD_INFO1] = p [0] ; 398 COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; 399 return (false) ; 400 } 401 402 /* === If no knobs, set default knobs =================================== */ 403 404 if (!knobs) 405 { 406 colamd_set_defaults (default_knobs) ; 407 knobs = default_knobs ; 408 } 409 410 /* === Allocate the Row and Col arrays from array A ===================== */ 411 412 Col_size = colamd_c (n_col) ; 413 Row_size = colamd_r (n_row) ; 414 need = 2*nnz + n_col + Col_size + Row_size ; 415 416 if (need > Alen) 417 { 418 /* not enough space in array A to perform the ordering */ 419 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; 420 stats [COLAMD_INFO1] = need ; 421 stats [COLAMD_INFO2] = Alen ; 422 COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); 423 return (false) ; 424 } 425 426 Alen -= Col_size + Row_size ; 427 Col = (colamd_col<IndexType> *) &A [Alen] ; 428 Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ; 429 430 /* === Construct the row and column data structures ===================== */ 431 432 if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) 433 { 434 /* input matrix is invalid */ 435 COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ; 436 return (false) ; 437 } 438 439 /* === Initialize scores, kill dense rows/columns ======================= */ 440 441 Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs, 442 &n_row2, &n_col2, &max_deg) ; 443 444 /* === Order the supercolumns =========================================== */ 445 446 ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p, 447 n_col2, max_deg, 2*nnz) ; 448 449 /* === Order the non-principal columns ================================== */ 450 451 Eigen::internal::order_children (n_col, Col, p) ; 452 453 /* === Return statistics in stats ======================================= */ 454 455 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; 456 stats [COLAMD_DENSE_COL] = n_col - n_col2 ; 457 stats [COLAMD_DEFRAG_COUNT] = ngarbage ; 458 COLAMD_DEBUG0 (("colamd: done.\n")) ; 459 return (true) ; 460} 461 462/* ========================================================================== */ 463/* === NON-USER-CALLABLE ROUTINES: ========================================== */ 464/* ========================================================================== */ 465 466/* There are no user-callable routines beyond this point in the file */ 467 468 469/* ========================================================================== */ 470/* === init_rows_cols ======================================================= */ 471/* ========================================================================== */ 472 473/* 474 Takes the column form of the matrix in A and creates the row form of the 475 matrix. Also, row and column attributes are stored in the Col and Row 476 structs. If the columns are un-sorted or contain duplicate row indices, 477 this routine will also sort and remove duplicate row indices from the 478 column form of the matrix. Returns false if the matrix is invalid, 479 true otherwise. Not user-callable. 480*/ 481template <typename IndexType> 482static IndexType init_rows_cols /* returns true if OK, or false otherwise */ 483 ( 484 /* === Parameters ======================================================= */ 485 486 IndexType n_row, /* number of rows of A */ 487 IndexType n_col, /* number of columns of A */ 488 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 489 colamd_col<IndexType> Col [], /* of size n_col+1 */ 490 IndexType A [], /* row indices of A, of size Alen */ 491 IndexType p [], /* pointers to columns in A, of size n_col+1 */ 492 IndexType stats [COLAMD_STATS] /* colamd statistics */ 493 ) 494{ 495 /* === Local variables ================================================== */ 496 497 IndexType col ; /* a column index */ 498 IndexType row ; /* a row index */ 499 IndexType *cp ; /* a column pointer */ 500 IndexType *cp_end ; /* a pointer to the end of a column */ 501 IndexType *rp ; /* a row pointer */ 502 IndexType *rp_end ; /* a pointer to the end of a row */ 503 IndexType last_row ; /* previous row */ 504 505 /* === Initialize columns, and check column pointers ==================== */ 506 507 for (col = 0 ; col < n_col ; col++) 508 { 509 Col [col].start = p [col] ; 510 Col [col].length = p [col+1] - p [col] ; 511 512 if ((Col [col].length) < 0) // extra parentheses to work-around gcc bug 10200 513 { 514 /* column pointers must be non-decreasing */ 515 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; 516 stats [COLAMD_INFO1] = col ; 517 stats [COLAMD_INFO2] = Col [col].length ; 518 COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; 519 return (false) ; 520 } 521 522 Col [col].shared1.thickness = 1 ; 523 Col [col].shared2.score = 0 ; 524 Col [col].shared3.prev = COLAMD_EMPTY ; 525 Col [col].shared4.degree_next = COLAMD_EMPTY ; 526 } 527 528 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ 529 530 /* === Scan columns, compute row degrees, and check row indices ========= */ 531 532 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ 533 534 for (row = 0 ; row < n_row ; row++) 535 { 536 Row [row].length = 0 ; 537 Row [row].shared2.mark = -1 ; 538 } 539 540 for (col = 0 ; col < n_col ; col++) 541 { 542 last_row = -1 ; 543 544 cp = &A [p [col]] ; 545 cp_end = &A [p [col+1]] ; 546 547 while (cp < cp_end) 548 { 549 row = *cp++ ; 550 551 /* make sure row indices within range */ 552 if (row < 0 || row >= n_row) 553 { 554 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; 555 stats [COLAMD_INFO1] = col ; 556 stats [COLAMD_INFO2] = row ; 557 stats [COLAMD_INFO3] = n_row ; 558 COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; 559 return (false) ; 560 } 561 562 if (row <= last_row || Row [row].shared2.mark == col) 563 { 564 /* row index are unsorted or repeated (or both), thus col */ 565 /* is jumbled. This is a notice, not an error condition. */ 566 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; 567 stats [COLAMD_INFO1] = col ; 568 stats [COLAMD_INFO2] = row ; 569 (stats [COLAMD_INFO3]) ++ ; 570 COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); 571 } 572 573 if (Row [row].shared2.mark != col) 574 { 575 Row [row].length++ ; 576 } 577 else 578 { 579 /* this is a repeated entry in the column, */ 580 /* it will be removed */ 581 Col [col].length-- ; 582 } 583 584 /* mark the row as having been seen in this column */ 585 Row [row].shared2.mark = col ; 586 587 last_row = row ; 588 } 589 } 590 591 /* === Compute row pointers ============================================= */ 592 593 /* row form of the matrix starts directly after the column */ 594 /* form of matrix in A */ 595 Row [0].start = p [n_col] ; 596 Row [0].shared1.p = Row [0].start ; 597 Row [0].shared2.mark = -1 ; 598 for (row = 1 ; row < n_row ; row++) 599 { 600 Row [row].start = Row [row-1].start + Row [row-1].length ; 601 Row [row].shared1.p = Row [row].start ; 602 Row [row].shared2.mark = -1 ; 603 } 604 605 /* === Create row form ================================================== */ 606 607 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 608 { 609 /* if cols jumbled, watch for repeated row indices */ 610 for (col = 0 ; col < n_col ; col++) 611 { 612 cp = &A [p [col]] ; 613 cp_end = &A [p [col+1]] ; 614 while (cp < cp_end) 615 { 616 row = *cp++ ; 617 if (Row [row].shared2.mark != col) 618 { 619 A [(Row [row].shared1.p)++] = col ; 620 Row [row].shared2.mark = col ; 621 } 622 } 623 } 624 } 625 else 626 { 627 /* if cols not jumbled, we don't need the mark (this is faster) */ 628 for (col = 0 ; col < n_col ; col++) 629 { 630 cp = &A [p [col]] ; 631 cp_end = &A [p [col+1]] ; 632 while (cp < cp_end) 633 { 634 A [(Row [*cp++].shared1.p)++] = col ; 635 } 636 } 637 } 638 639 /* === Clear the row marks and set row degrees ========================== */ 640 641 for (row = 0 ; row < n_row ; row++) 642 { 643 Row [row].shared2.mark = 0 ; 644 Row [row].shared1.degree = Row [row].length ; 645 } 646 647 /* === See if we need to re-create columns ============================== */ 648 649 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 650 { 651 COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; 652 653 654 /* === Compute col pointers ========================================= */ 655 656 /* col form of the matrix starts at A [0]. */ 657 /* Note, we may have a gap between the col form and the row */ 658 /* form if there were duplicate entries, if so, it will be */ 659 /* removed upon the first garbage collection */ 660 Col [0].start = 0 ; 661 p [0] = Col [0].start ; 662 for (col = 1 ; col < n_col ; col++) 663 { 664 /* note that the lengths here are for pruned columns, i.e. */ 665 /* no duplicate row indices will exist for these columns */ 666 Col [col].start = Col [col-1].start + Col [col-1].length ; 667 p [col] = Col [col].start ; 668 } 669 670 /* === Re-create col form =========================================== */ 671 672 for (row = 0 ; row < n_row ; row++) 673 { 674 rp = &A [Row [row].start] ; 675 rp_end = rp + Row [row].length ; 676 while (rp < rp_end) 677 { 678 A [(p [*rp++])++] = row ; 679 } 680 } 681 } 682 683 /* === Done. Matrix is not (or no longer) jumbled ====================== */ 684 685 return (true) ; 686} 687 688 689/* ========================================================================== */ 690/* === init_scoring ========================================================= */ 691/* ========================================================================== */ 692 693/* 694 Kills dense or empty columns and rows, calculates an initial score for 695 each column, and places all columns in the degree lists. Not user-callable. 696*/ 697template <typename IndexType> 698static void init_scoring 699 ( 700 /* === Parameters ======================================================= */ 701 702 IndexType n_row, /* number of rows of A */ 703 IndexType n_col, /* number of columns of A */ 704 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 705 colamd_col<IndexType> Col [], /* of size n_col+1 */ 706 IndexType A [], /* column form and row form of A */ 707 IndexType head [], /* of size n_col+1 */ 708 double knobs [COLAMD_KNOBS],/* parameters */ 709 IndexType *p_n_row2, /* number of non-dense, non-empty rows */ 710 IndexType *p_n_col2, /* number of non-dense, non-empty columns */ 711 IndexType *p_max_deg /* maximum row degree */ 712 ) 713{ 714 /* === Local variables ================================================== */ 715 716 IndexType c ; /* a column index */ 717 IndexType r, row ; /* a row index */ 718 IndexType *cp ; /* a column pointer */ 719 IndexType deg ; /* degree of a row or column */ 720 IndexType *cp_end ; /* a pointer to the end of a column */ 721 IndexType *new_cp ; /* new column pointer */ 722 IndexType col_length ; /* length of pruned column */ 723 IndexType score ; /* current column score */ 724 IndexType n_col2 ; /* number of non-dense, non-empty columns */ 725 IndexType n_row2 ; /* number of non-dense, non-empty rows */ 726 IndexType dense_row_count ; /* remove rows with more entries than this */ 727 IndexType dense_col_count ; /* remove cols with more entries than this */ 728 IndexType min_score ; /* smallest column score */ 729 IndexType max_deg ; /* maximum row degree */ 730 IndexType next_col ; /* Used to add to degree list.*/ 731 732 733 /* === Extract knobs ==================================================== */ 734 735 dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ; 736 dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ; 737 COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; 738 max_deg = 0 ; 739 n_col2 = n_col ; 740 n_row2 = n_row ; 741 742 /* === Kill empty columns =============================================== */ 743 744 /* Put the empty columns at the end in their natural order, so that LU */ 745 /* factorization can proceed as far as possible. */ 746 for (c = n_col-1 ; c >= 0 ; c--) 747 { 748 deg = Col [c].length ; 749 if (deg == 0) 750 { 751 /* this is a empty column, kill and order it last */ 752 Col [c].shared2.order = --n_col2 ; 753 KILL_PRINCIPAL_COL (c) ; 754 } 755 } 756 COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; 757 758 /* === Kill dense columns =============================================== */ 759 760 /* Put the dense columns at the end, in their natural order */ 761 for (c = n_col-1 ; c >= 0 ; c--) 762 { 763 /* skip any dead columns */ 764 if (COL_IS_DEAD (c)) 765 { 766 continue ; 767 } 768 deg = Col [c].length ; 769 if (deg > dense_col_count) 770 { 771 /* this is a dense column, kill and order it last */ 772 Col [c].shared2.order = --n_col2 ; 773 /* decrement the row degrees */ 774 cp = &A [Col [c].start] ; 775 cp_end = cp + Col [c].length ; 776 while (cp < cp_end) 777 { 778 Row [*cp++].shared1.degree-- ; 779 } 780 KILL_PRINCIPAL_COL (c) ; 781 } 782 } 783 COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; 784 785 /* === Kill dense and empty rows ======================================== */ 786 787 for (r = 0 ; r < n_row ; r++) 788 { 789 deg = Row [r].shared1.degree ; 790 COLAMD_ASSERT (deg >= 0 && deg <= n_col) ; 791 if (deg > dense_row_count || deg == 0) 792 { 793 /* kill a dense or empty row */ 794 KILL_ROW (r) ; 795 --n_row2 ; 796 } 797 else 798 { 799 /* keep track of max degree of remaining rows */ 800 max_deg = numext::maxi(max_deg, deg) ; 801 } 802 } 803 COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; 804 805 /* === Compute initial column scores ==================================== */ 806 807 /* At this point the row degrees are accurate. They reflect the number */ 808 /* of "live" (non-dense) columns in each row. No empty rows exist. */ 809 /* Some "live" columns may contain only dead rows, however. These are */ 810 /* pruned in the code below. */ 811 812 /* now find the initial matlab score for each column */ 813 for (c = n_col-1 ; c >= 0 ; c--) 814 { 815 /* skip dead column */ 816 if (COL_IS_DEAD (c)) 817 { 818 continue ; 819 } 820 score = 0 ; 821 cp = &A [Col [c].start] ; 822 new_cp = cp ; 823 cp_end = cp + Col [c].length ; 824 while (cp < cp_end) 825 { 826 /* get a row */ 827 row = *cp++ ; 828 /* skip if dead */ 829 if (ROW_IS_DEAD (row)) 830 { 831 continue ; 832 } 833 /* compact the column */ 834 *new_cp++ = row ; 835 /* add row's external degree */ 836 score += Row [row].shared1.degree - 1 ; 837 /* guard against integer overflow */ 838 score = numext::mini(score, n_col) ; 839 } 840 /* determine pruned column length */ 841 col_length = (IndexType) (new_cp - &A [Col [c].start]) ; 842 if (col_length == 0) 843 { 844 /* a newly-made null column (all rows in this col are "dense" */ 845 /* and have already been killed) */ 846 COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ; 847 Col [c].shared2.order = --n_col2 ; 848 KILL_PRINCIPAL_COL (c) ; 849 } 850 else 851 { 852 /* set column length and set score */ 853 COLAMD_ASSERT (score >= 0) ; 854 COLAMD_ASSERT (score <= n_col) ; 855 Col [c].length = col_length ; 856 Col [c].shared2.score = score ; 857 } 858 } 859 COLAMD_DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", 860 n_col-n_col2)) ; 861 862 /* At this point, all empty rows and columns are dead. All live columns */ 863 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ 864 /* yet). Rows may contain dead columns, but all live rows contain at */ 865 /* least one live column. */ 866 867 /* === Initialize degree lists ========================================== */ 868 869 870 /* clear the hash buckets */ 871 for (c = 0 ; c <= n_col ; c++) 872 { 873 head [c] = COLAMD_EMPTY ; 874 } 875 min_score = n_col ; 876 /* place in reverse order, so low column indices are at the front */ 877 /* of the lists. This is to encourage natural tie-breaking */ 878 for (c = n_col-1 ; c >= 0 ; c--) 879 { 880 /* only add principal columns to degree lists */ 881 if (COL_IS_ALIVE (c)) 882 { 883 COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n", 884 c, Col [c].shared2.score, min_score, n_col)) ; 885 886 /* === Add columns score to DList =============================== */ 887 888 score = Col [c].shared2.score ; 889 890 COLAMD_ASSERT (min_score >= 0) ; 891 COLAMD_ASSERT (min_score <= n_col) ; 892 COLAMD_ASSERT (score >= 0) ; 893 COLAMD_ASSERT (score <= n_col) ; 894 COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ; 895 896 /* now add this column to dList at proper score location */ 897 next_col = head [score] ; 898 Col [c].shared3.prev = COLAMD_EMPTY ; 899 Col [c].shared4.degree_next = next_col ; 900 901 /* if there already was a column with the same score, set its */ 902 /* previous pointer to this new column */ 903 if (next_col != COLAMD_EMPTY) 904 { 905 Col [next_col].shared3.prev = c ; 906 } 907 head [score] = c ; 908 909 /* see if this score is less than current min */ 910 min_score = numext::mini(min_score, score) ; 911 912 913 } 914 } 915 916 917 /* === Return number of remaining columns, and max row degree =========== */ 918 919 *p_n_col2 = n_col2 ; 920 *p_n_row2 = n_row2 ; 921 *p_max_deg = max_deg ; 922} 923 924 925/* ========================================================================== */ 926/* === find_ordering ======================================================== */ 927/* ========================================================================== */ 928 929/* 930 Order the principal columns of the supercolumn form of the matrix 931 (no supercolumns on input). Uses a minimum approximate column minimum 932 degree ordering method. Not user-callable. 933*/ 934template <typename IndexType> 935static IndexType find_ordering /* return the number of garbage collections */ 936 ( 937 /* === Parameters ======================================================= */ 938 939 IndexType n_row, /* number of rows of A */ 940 IndexType n_col, /* number of columns of A */ 941 IndexType Alen, /* size of A, 2*nnz + n_col or larger */ 942 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 943 colamd_col<IndexType> Col [], /* of size n_col+1 */ 944 IndexType A [], /* column form and row form of A */ 945 IndexType head [], /* of size n_col+1 */ 946 IndexType n_col2, /* Remaining columns to order */ 947 IndexType max_deg, /* Maximum row degree */ 948 IndexType pfree /* index of first free slot (2*nnz on entry) */ 949 ) 950{ 951 /* === Local variables ================================================== */ 952 953 IndexType k ; /* current pivot ordering step */ 954 IndexType pivot_col ; /* current pivot column */ 955 IndexType *cp ; /* a column pointer */ 956 IndexType *rp ; /* a row pointer */ 957 IndexType pivot_row ; /* current pivot row */ 958 IndexType *new_cp ; /* modified column pointer */ 959 IndexType *new_rp ; /* modified row pointer */ 960 IndexType pivot_row_start ; /* pointer to start of pivot row */ 961 IndexType pivot_row_degree ; /* number of columns in pivot row */ 962 IndexType pivot_row_length ; /* number of supercolumns in pivot row */ 963 IndexType pivot_col_score ; /* score of pivot column */ 964 IndexType needed_memory ; /* free space needed for pivot row */ 965 IndexType *cp_end ; /* pointer to the end of a column */ 966 IndexType *rp_end ; /* pointer to the end of a row */ 967 IndexType row ; /* a row index */ 968 IndexType col ; /* a column index */ 969 IndexType max_score ; /* maximum possible score */ 970 IndexType cur_score ; /* score of current column */ 971 unsigned int hash ; /* hash value for supernode detection */ 972 IndexType head_column ; /* head of hash bucket */ 973 IndexType first_col ; /* first column in hash bucket */ 974 IndexType tag_mark ; /* marker value for mark array */ 975 IndexType row_mark ; /* Row [row].shared2.mark */ 976 IndexType set_difference ; /* set difference size of row with pivot row */ 977 IndexType min_score ; /* smallest column score */ 978 IndexType col_thickness ; /* "thickness" (no. of columns in a supercol) */ 979 IndexType max_mark ; /* maximum value of tag_mark */ 980 IndexType pivot_col_thickness ; /* number of columns represented by pivot col */ 981 IndexType prev_col ; /* Used by Dlist operations. */ 982 IndexType next_col ; /* Used by Dlist operations. */ 983 IndexType ngarbage ; /* number of garbage collections performed */ 984 985 986 /* === Initialization and clear mark ==================================== */ 987 988 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */ 989 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 990 min_score = 0 ; 991 ngarbage = 0 ; 992 COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; 993 994 /* === Order the columns ================================================ */ 995 996 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) 997 { 998 999 /* === Select pivot column, and order it ============================ */ 1000 1001 /* make sure degree list isn't empty */ 1002 COLAMD_ASSERT (min_score >= 0) ; 1003 COLAMD_ASSERT (min_score <= n_col) ; 1004 COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ; 1005 1006 /* get pivot column from head of minimum degree list */ 1007 while (min_score < n_col && head [min_score] == COLAMD_EMPTY) 1008 { 1009 min_score++ ; 1010 } 1011 pivot_col = head [min_score] ; 1012 COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; 1013 next_col = Col [pivot_col].shared4.degree_next ; 1014 head [min_score] = next_col ; 1015 if (next_col != COLAMD_EMPTY) 1016 { 1017 Col [next_col].shared3.prev = COLAMD_EMPTY ; 1018 } 1019 1020 COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ; 1021 COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ; 1022 1023 /* remember score for defrag check */ 1024 pivot_col_score = Col [pivot_col].shared2.score ; 1025 1026 /* the pivot column is the kth column in the pivot order */ 1027 Col [pivot_col].shared2.order = k ; 1028 1029 /* increment order count by column thickness */ 1030 pivot_col_thickness = Col [pivot_col].shared1.thickness ; 1031 k += pivot_col_thickness ; 1032 COLAMD_ASSERT (pivot_col_thickness > 0) ; 1033 1034 /* === Garbage_collection, if necessary ============================= */ 1035 1036 needed_memory = numext::mini(pivot_col_score, n_col - k) ; 1037 if (pfree + needed_memory >= Alen) 1038 { 1039 pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; 1040 ngarbage++ ; 1041 /* after garbage collection we will have enough */ 1042 COLAMD_ASSERT (pfree + needed_memory < Alen) ; 1043 /* garbage collection has wiped out the Row[].shared2.mark array */ 1044 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 1045 1046 } 1047 1048 /* === Compute pivot row pattern ==================================== */ 1049 1050 /* get starting location for this new merged row */ 1051 pivot_row_start = pfree ; 1052 1053 /* initialize new row counts to zero */ 1054 pivot_row_degree = 0 ; 1055 1056 /* tag pivot column as having been visited so it isn't included */ 1057 /* in merged pivot row */ 1058 Col [pivot_col].shared1.thickness = -pivot_col_thickness ; 1059 1060 /* pivot row is the union of all rows in the pivot column pattern */ 1061 cp = &A [Col [pivot_col].start] ; 1062 cp_end = cp + Col [pivot_col].length ; 1063 while (cp < cp_end) 1064 { 1065 /* get a row */ 1066 row = *cp++ ; 1067 COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; 1068 /* skip if row is dead */ 1069 if (ROW_IS_DEAD (row)) 1070 { 1071 continue ; 1072 } 1073 rp = &A [Row [row].start] ; 1074 rp_end = rp + Row [row].length ; 1075 while (rp < rp_end) 1076 { 1077 /* get a column */ 1078 col = *rp++ ; 1079 /* add the column, if alive and untagged */ 1080 col_thickness = Col [col].shared1.thickness ; 1081 if (col_thickness > 0 && COL_IS_ALIVE (col)) 1082 { 1083 /* tag column in pivot row */ 1084 Col [col].shared1.thickness = -col_thickness ; 1085 COLAMD_ASSERT (pfree < Alen) ; 1086 /* place column in pivot row */ 1087 A [pfree++] = col ; 1088 pivot_row_degree += col_thickness ; 1089 } 1090 } 1091 } 1092 1093 /* clear tag on pivot column */ 1094 Col [pivot_col].shared1.thickness = pivot_col_thickness ; 1095 max_deg = numext::maxi(max_deg, pivot_row_degree) ; 1096 1097 1098 /* === Kill all rows used to construct pivot row ==================== */ 1099 1100 /* also kill pivot row, temporarily */ 1101 cp = &A [Col [pivot_col].start] ; 1102 cp_end = cp + Col [pivot_col].length ; 1103 while (cp < cp_end) 1104 { 1105 /* may be killing an already dead row */ 1106 row = *cp++ ; 1107 COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ; 1108 KILL_ROW (row) ; 1109 } 1110 1111 /* === Select a row index to use as the new pivot row =============== */ 1112 1113 pivot_row_length = pfree - pivot_row_start ; 1114 if (pivot_row_length > 0) 1115 { 1116 /* pick the "pivot" row arbitrarily (first row in col) */ 1117 pivot_row = A [Col [pivot_col].start] ; 1118 COLAMD_DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; 1119 } 1120 else 1121 { 1122 /* there is no pivot row, since it is of zero length */ 1123 pivot_row = COLAMD_EMPTY ; 1124 COLAMD_ASSERT (pivot_row_length == 0) ; 1125 } 1126 COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; 1127 1128 /* === Approximate degree computation =============================== */ 1129 1130 /* Here begins the computation of the approximate degree. The column */ 1131 /* score is the sum of the pivot row "length", plus the size of the */ 1132 /* set differences of each row in the column minus the pattern of the */ 1133 /* pivot row itself. The column ("thickness") itself is also */ 1134 /* excluded from the column score (we thus use an approximate */ 1135 /* external degree). */ 1136 1137 /* The time taken by the following code (compute set differences, and */ 1138 /* add them up) is proportional to the size of the data structure */ 1139 /* being scanned - that is, the sum of the sizes of each column in */ 1140 /* the pivot row. Thus, the amortized time to compute a column score */ 1141 /* is proportional to the size of that column (where size, in this */ 1142 /* context, is the column "length", or the number of row indices */ 1143 /* in that column). The number of row indices in a column is */ 1144 /* monotonically non-decreasing, from the length of the original */ 1145 /* column on input to colamd. */ 1146 1147 /* === Compute set differences ====================================== */ 1148 1149 COLAMD_DEBUG3 (("** Computing set differences phase. **\n")) ; 1150 1151 /* pivot row is currently dead - it will be revived later. */ 1152 1153 COLAMD_DEBUG3 (("Pivot row: ")) ; 1154 /* for each column in pivot row */ 1155 rp = &A [pivot_row_start] ; 1156 rp_end = rp + pivot_row_length ; 1157 while (rp < rp_end) 1158 { 1159 col = *rp++ ; 1160 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 1161 COLAMD_DEBUG3 (("Col: %d\n", col)) ; 1162 1163 /* clear tags used to construct pivot row pattern */ 1164 col_thickness = -Col [col].shared1.thickness ; 1165 COLAMD_ASSERT (col_thickness > 0) ; 1166 Col [col].shared1.thickness = col_thickness ; 1167 1168 /* === Remove column from degree list =========================== */ 1169 1170 cur_score = Col [col].shared2.score ; 1171 prev_col = Col [col].shared3.prev ; 1172 next_col = Col [col].shared4.degree_next ; 1173 COLAMD_ASSERT (cur_score >= 0) ; 1174 COLAMD_ASSERT (cur_score <= n_col) ; 1175 COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ; 1176 if (prev_col == COLAMD_EMPTY) 1177 { 1178 head [cur_score] = next_col ; 1179 } 1180 else 1181 { 1182 Col [prev_col].shared4.degree_next = next_col ; 1183 } 1184 if (next_col != COLAMD_EMPTY) 1185 { 1186 Col [next_col].shared3.prev = prev_col ; 1187 } 1188 1189 /* === Scan the column ========================================== */ 1190 1191 cp = &A [Col [col].start] ; 1192 cp_end = cp + Col [col].length ; 1193 while (cp < cp_end) 1194 { 1195 /* get a row */ 1196 row = *cp++ ; 1197 row_mark = Row [row].shared2.mark ; 1198 /* skip if dead */ 1199 if (ROW_IS_MARKED_DEAD (row_mark)) 1200 { 1201 continue ; 1202 } 1203 COLAMD_ASSERT (row != pivot_row) ; 1204 set_difference = row_mark - tag_mark ; 1205 /* check if the row has been seen yet */ 1206 if (set_difference < 0) 1207 { 1208 COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ; 1209 set_difference = Row [row].shared1.degree ; 1210 } 1211 /* subtract column thickness from this row's set difference */ 1212 set_difference -= col_thickness ; 1213 COLAMD_ASSERT (set_difference >= 0) ; 1214 /* absorb this row if the set difference becomes zero */ 1215 if (set_difference == 0) 1216 { 1217 COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; 1218 KILL_ROW (row) ; 1219 } 1220 else 1221 { 1222 /* save the new mark */ 1223 Row [row].shared2.mark = set_difference + tag_mark ; 1224 } 1225 } 1226 } 1227 1228 1229 /* === Add up set differences for each column ======================= */ 1230 1231 COLAMD_DEBUG3 (("** Adding set differences phase. **\n")) ; 1232 1233 /* for each column in pivot row */ 1234 rp = &A [pivot_row_start] ; 1235 rp_end = rp + pivot_row_length ; 1236 while (rp < rp_end) 1237 { 1238 /* get a column */ 1239 col = *rp++ ; 1240 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 1241 hash = 0 ; 1242 cur_score = 0 ; 1243 cp = &A [Col [col].start] ; 1244 /* compact the column */ 1245 new_cp = cp ; 1246 cp_end = cp + Col [col].length ; 1247 1248 COLAMD_DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; 1249 1250 while (cp < cp_end) 1251 { 1252 /* get a row */ 1253 row = *cp++ ; 1254 COLAMD_ASSERT(row >= 0 && row < n_row) ; 1255 row_mark = Row [row].shared2.mark ; 1256 /* skip if dead */ 1257 if (ROW_IS_MARKED_DEAD (row_mark)) 1258 { 1259 continue ; 1260 } 1261 COLAMD_ASSERT (row_mark > tag_mark) ; 1262 /* compact the column */ 1263 *new_cp++ = row ; 1264 /* compute hash function */ 1265 hash += row ; 1266 /* add set difference */ 1267 cur_score += row_mark - tag_mark ; 1268 /* integer overflow... */ 1269 cur_score = numext::mini(cur_score, n_col) ; 1270 } 1271 1272 /* recompute the column's length */ 1273 Col [col].length = (IndexType) (new_cp - &A [Col [col].start]) ; 1274 1275 /* === Further mass elimination ================================= */ 1276 1277 if (Col [col].length == 0) 1278 { 1279 COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ; 1280 /* nothing left but the pivot row in this column */ 1281 KILL_PRINCIPAL_COL (col) ; 1282 pivot_row_degree -= Col [col].shared1.thickness ; 1283 COLAMD_ASSERT (pivot_row_degree >= 0) ; 1284 /* order it */ 1285 Col [col].shared2.order = k ; 1286 /* increment order count by column thickness */ 1287 k += Col [col].shared1.thickness ; 1288 } 1289 else 1290 { 1291 /* === Prepare for supercolumn detection ==================== */ 1292 1293 COLAMD_DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; 1294 1295 /* save score so far */ 1296 Col [col].shared2.score = cur_score ; 1297 1298 /* add column to hash table, for supercolumn detection */ 1299 hash %= n_col + 1 ; 1300 1301 COLAMD_DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; 1302 COLAMD_ASSERT (hash <= n_col) ; 1303 1304 head_column = head [hash] ; 1305 if (head_column > COLAMD_EMPTY) 1306 { 1307 /* degree list "hash" is non-empty, use prev (shared3) of */ 1308 /* first column in degree list as head of hash bucket */ 1309 first_col = Col [head_column].shared3.headhash ; 1310 Col [head_column].shared3.headhash = col ; 1311 } 1312 else 1313 { 1314 /* degree list "hash" is empty, use head as hash bucket */ 1315 first_col = - (head_column + 2) ; 1316 head [hash] = - (col + 2) ; 1317 } 1318 Col [col].shared4.hash_next = first_col ; 1319 1320 /* save hash function in Col [col].shared3.hash */ 1321 Col [col].shared3.hash = (IndexType) hash ; 1322 COLAMD_ASSERT (COL_IS_ALIVE (col)) ; 1323 } 1324 } 1325 1326 /* The approximate external column degree is now computed. */ 1327 1328 /* === Supercolumn detection ======================================== */ 1329 1330 COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ; 1331 1332 Eigen::internal::detect_super_cols (Col, A, head, pivot_row_start, pivot_row_length) ; 1333 1334 /* === Kill the pivotal column ====================================== */ 1335 1336 KILL_PRINCIPAL_COL (pivot_col) ; 1337 1338 /* === Clear mark =================================================== */ 1339 1340 tag_mark += (max_deg + 1) ; 1341 if (tag_mark >= max_mark) 1342 { 1343 COLAMD_DEBUG2 (("clearing tag_mark\n")) ; 1344 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 1345 } 1346 1347 /* === Finalize the new pivot row, and column scores ================ */ 1348 1349 COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ; 1350 1351 /* for each column in pivot row */ 1352 rp = &A [pivot_row_start] ; 1353 /* compact the pivot row */ 1354 new_rp = rp ; 1355 rp_end = rp + pivot_row_length ; 1356 while (rp < rp_end) 1357 { 1358 col = *rp++ ; 1359 /* skip dead columns */ 1360 if (COL_IS_DEAD (col)) 1361 { 1362 continue ; 1363 } 1364 *new_rp++ = col ; 1365 /* add new pivot row to column */ 1366 A [Col [col].start + (Col [col].length++)] = pivot_row ; 1367 1368 /* retrieve score so far and add on pivot row's degree. */ 1369 /* (we wait until here for this in case the pivot */ 1370 /* row's degree was reduced due to mass elimination). */ 1371 cur_score = Col [col].shared2.score + pivot_row_degree ; 1372 1373 /* calculate the max possible score as the number of */ 1374 /* external columns minus the 'k' value minus the */ 1375 /* columns thickness */ 1376 max_score = n_col - k - Col [col].shared1.thickness ; 1377 1378 /* make the score the external degree of the union-of-rows */ 1379 cur_score -= Col [col].shared1.thickness ; 1380 1381 /* make sure score is less or equal than the max score */ 1382 cur_score = numext::mini(cur_score, max_score) ; 1383 COLAMD_ASSERT (cur_score >= 0) ; 1384 1385 /* store updated score */ 1386 Col [col].shared2.score = cur_score ; 1387 1388 /* === Place column back in degree list ========================= */ 1389 1390 COLAMD_ASSERT (min_score >= 0) ; 1391 COLAMD_ASSERT (min_score <= n_col) ; 1392 COLAMD_ASSERT (cur_score >= 0) ; 1393 COLAMD_ASSERT (cur_score <= n_col) ; 1394 COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ; 1395 next_col = head [cur_score] ; 1396 Col [col].shared4.degree_next = next_col ; 1397 Col [col].shared3.prev = COLAMD_EMPTY ; 1398 if (next_col != COLAMD_EMPTY) 1399 { 1400 Col [next_col].shared3.prev = col ; 1401 } 1402 head [cur_score] = col ; 1403 1404 /* see if this score is less than current min */ 1405 min_score = numext::mini(min_score, cur_score) ; 1406 1407 } 1408 1409 /* === Resurrect the new pivot row ================================== */ 1410 1411 if (pivot_row_degree > 0) 1412 { 1413 /* update pivot row length to reflect any cols that were killed */ 1414 /* during super-col detection and mass elimination */ 1415 Row [pivot_row].start = pivot_row_start ; 1416 Row [pivot_row].length = (IndexType) (new_rp - &A[pivot_row_start]) ; 1417 Row [pivot_row].shared1.degree = pivot_row_degree ; 1418 Row [pivot_row].shared2.mark = 0 ; 1419 /* pivot row is no longer dead */ 1420 } 1421 } 1422 1423 /* === All principal columns have now been ordered ====================== */ 1424 1425 return (ngarbage) ; 1426} 1427 1428 1429/* ========================================================================== */ 1430/* === order_children ======================================================= */ 1431/* ========================================================================== */ 1432 1433/* 1434 The find_ordering routine has ordered all of the principal columns (the 1435 representatives of the supercolumns). The non-principal columns have not 1436 yet been ordered. This routine orders those columns by walking up the 1437 parent tree (a column is a child of the column which absorbed it). The 1438 final permutation vector is then placed in p [0 ... n_col-1], with p [0] 1439 being the first column, and p [n_col-1] being the last. It doesn't look 1440 like it at first glance, but be assured that this routine takes time linear 1441 in the number of columns. Although not immediately obvious, the time 1442 taken by this routine is O (n_col), that is, linear in the number of 1443 columns. Not user-callable. 1444*/ 1445template <typename IndexType> 1446static inline void order_children 1447( 1448 /* === Parameters ======================================================= */ 1449 1450 IndexType n_col, /* number of columns of A */ 1451 colamd_col<IndexType> Col [], /* of size n_col+1 */ 1452 IndexType p [] /* p [0 ... n_col-1] is the column permutation*/ 1453 ) 1454{ 1455 /* === Local variables ================================================== */ 1456 1457 IndexType i ; /* loop counter for all columns */ 1458 IndexType c ; /* column index */ 1459 IndexType parent ; /* index of column's parent */ 1460 IndexType order ; /* column's order */ 1461 1462 /* === Order each non-principal column ================================== */ 1463 1464 for (i = 0 ; i < n_col ; i++) 1465 { 1466 /* find an un-ordered non-principal column */ 1467 COLAMD_ASSERT (COL_IS_DEAD (i)) ; 1468 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY) 1469 { 1470 parent = i ; 1471 /* once found, find its principal parent */ 1472 do 1473 { 1474 parent = Col [parent].shared1.parent ; 1475 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; 1476 1477 /* now, order all un-ordered non-principal columns along path */ 1478 /* to this parent. collapse tree at the same time */ 1479 c = i ; 1480 /* get order of parent */ 1481 order = Col [parent].shared2.order ; 1482 1483 do 1484 { 1485 COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ; 1486 1487 /* order this column */ 1488 Col [c].shared2.order = order++ ; 1489 /* collaps tree */ 1490 Col [c].shared1.parent = parent ; 1491 1492 /* get immediate parent of this column */ 1493 c = Col [c].shared1.parent ; 1494 1495 /* continue until we hit an ordered column. There are */ 1496 /* guarranteed not to be anymore unordered columns */ 1497 /* above an ordered column */ 1498 } while (Col [c].shared2.order == COLAMD_EMPTY) ; 1499 1500 /* re-order the super_col parent to largest order for this group */ 1501 Col [parent].shared2.order = order ; 1502 } 1503 } 1504 1505 /* === Generate the permutation ========================================= */ 1506 1507 for (c = 0 ; c < n_col ; c++) 1508 { 1509 p [Col [c].shared2.order] = c ; 1510 } 1511} 1512 1513 1514/* ========================================================================== */ 1515/* === detect_super_cols ==================================================== */ 1516/* ========================================================================== */ 1517 1518/* 1519 Detects supercolumns by finding matches between columns in the hash buckets. 1520 Check amongst columns in the set A [row_start ... row_start + row_length-1]. 1521 The columns under consideration are currently *not* in the degree lists, 1522 and have already been placed in the hash buckets. 1523 1524 The hash bucket for columns whose hash function is equal to h is stored 1525 as follows: 1526 1527 if head [h] is >= 0, then head [h] contains a degree list, so: 1528 1529 head [h] is the first column in degree bucket h. 1530 Col [head [h]].headhash gives the first column in hash bucket h. 1531 1532 otherwise, the degree list is empty, and: 1533 1534 -(head [h] + 2) is the first column in hash bucket h. 1535 1536 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous 1537 column" pointer. Col [c].shared3.hash is used instead as the hash number 1538 for that column. The value of Col [c].shared4.hash_next is the next column 1539 in the same hash bucket. 1540 1541 Assuming no, or "few" hash collisions, the time taken by this routine is 1542 linear in the sum of the sizes (lengths) of each column whose score has 1543 just been computed in the approximate degree computation. 1544 Not user-callable. 1545*/ 1546template <typename IndexType> 1547static void detect_super_cols 1548( 1549 /* === Parameters ======================================================= */ 1550 1551 colamd_col<IndexType> Col [], /* of size n_col+1 */ 1552 IndexType A [], /* row indices of A */ 1553 IndexType head [], /* head of degree lists and hash buckets */ 1554 IndexType row_start, /* pointer to set of columns to check */ 1555 IndexType row_length /* number of columns to check */ 1556) 1557{ 1558 /* === Local variables ================================================== */ 1559 1560 IndexType hash ; /* hash value for a column */ 1561 IndexType *rp ; /* pointer to a row */ 1562 IndexType c ; /* a column index */ 1563 IndexType super_c ; /* column index of the column to absorb into */ 1564 IndexType *cp1 ; /* column pointer for column super_c */ 1565 IndexType *cp2 ; /* column pointer for column c */ 1566 IndexType length ; /* length of column super_c */ 1567 IndexType prev_c ; /* column preceding c in hash bucket */ 1568 IndexType i ; /* loop counter */ 1569 IndexType *rp_end ; /* pointer to the end of the row */ 1570 IndexType col ; /* a column index in the row to check */ 1571 IndexType head_column ; /* first column in hash bucket or degree list */ 1572 IndexType first_col ; /* first column in hash bucket */ 1573 1574 /* === Consider each column in the row ================================== */ 1575 1576 rp = &A [row_start] ; 1577 rp_end = rp + row_length ; 1578 while (rp < rp_end) 1579 { 1580 col = *rp++ ; 1581 if (COL_IS_DEAD (col)) 1582 { 1583 continue ; 1584 } 1585 1586 /* get hash number for this column */ 1587 hash = Col [col].shared3.hash ; 1588 COLAMD_ASSERT (hash <= n_col) ; 1589 1590 /* === Get the first column in this hash bucket ===================== */ 1591 1592 head_column = head [hash] ; 1593 if (head_column > COLAMD_EMPTY) 1594 { 1595 first_col = Col [head_column].shared3.headhash ; 1596 } 1597 else 1598 { 1599 first_col = - (head_column + 2) ; 1600 } 1601 1602 /* === Consider each column in the hash bucket ====================== */ 1603 1604 for (super_c = first_col ; super_c != COLAMD_EMPTY ; 1605 super_c = Col [super_c].shared4.hash_next) 1606 { 1607 COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ; 1608 COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ; 1609 length = Col [super_c].length ; 1610 1611 /* prev_c is the column preceding column c in the hash bucket */ 1612 prev_c = super_c ; 1613 1614 /* === Compare super_c with all columns after it ================ */ 1615 1616 for (c = Col [super_c].shared4.hash_next ; 1617 c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next) 1618 { 1619 COLAMD_ASSERT (c != super_c) ; 1620 COLAMD_ASSERT (COL_IS_ALIVE (c)) ; 1621 COLAMD_ASSERT (Col [c].shared3.hash == hash) ; 1622 1623 /* not identical if lengths or scores are different */ 1624 if (Col [c].length != length || 1625 Col [c].shared2.score != Col [super_c].shared2.score) 1626 { 1627 prev_c = c ; 1628 continue ; 1629 } 1630 1631 /* compare the two columns */ 1632 cp1 = &A [Col [super_c].start] ; 1633 cp2 = &A [Col [c].start] ; 1634 1635 for (i = 0 ; i < length ; i++) 1636 { 1637 /* the columns are "clean" (no dead rows) */ 1638 COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ; 1639 COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ; 1640 /* row indices will same order for both supercols, */ 1641 /* no gather scatter nessasary */ 1642 if (*cp1++ != *cp2++) 1643 { 1644 break ; 1645 } 1646 } 1647 1648 /* the two columns are different if the for-loop "broke" */ 1649 if (i != length) 1650 { 1651 prev_c = c ; 1652 continue ; 1653 } 1654 1655 /* === Got it! two columns are identical =================== */ 1656 1657 COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; 1658 1659 Col [super_c].shared1.thickness += Col [c].shared1.thickness ; 1660 Col [c].shared1.parent = super_c ; 1661 KILL_NON_PRINCIPAL_COL (c) ; 1662 /* order c later, in order_children() */ 1663 Col [c].shared2.order = COLAMD_EMPTY ; 1664 /* remove c from hash bucket */ 1665 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; 1666 } 1667 } 1668 1669 /* === Empty this hash bucket ======================================= */ 1670 1671 if (head_column > COLAMD_EMPTY) 1672 { 1673 /* corresponding degree list "hash" is not empty */ 1674 Col [head_column].shared3.headhash = COLAMD_EMPTY ; 1675 } 1676 else 1677 { 1678 /* corresponding degree list "hash" is empty */ 1679 head [hash] = COLAMD_EMPTY ; 1680 } 1681 } 1682} 1683 1684 1685/* ========================================================================== */ 1686/* === garbage_collection =================================================== */ 1687/* ========================================================================== */ 1688 1689/* 1690 Defragments and compacts columns and rows in the workspace A. Used when 1691 all avaliable memory has been used while performing row merging. Returns 1692 the index of the first free position in A, after garbage collection. The 1693 time taken by this routine is linear is the size of the array A, which is 1694 itself linear in the number of nonzeros in the input matrix. 1695 Not user-callable. 1696*/ 1697template <typename IndexType> 1698static IndexType garbage_collection /* returns the new value of pfree */ 1699 ( 1700 /* === Parameters ======================================================= */ 1701 1702 IndexType n_row, /* number of rows */ 1703 IndexType n_col, /* number of columns */ 1704 Colamd_Row<IndexType> Row [], /* row info */ 1705 colamd_col<IndexType> Col [], /* column info */ 1706 IndexType A [], /* A [0 ... Alen-1] holds the matrix */ 1707 IndexType *pfree /* &A [0] ... pfree is in use */ 1708 ) 1709{ 1710 /* === Local variables ================================================== */ 1711 1712 IndexType *psrc ; /* source pointer */ 1713 IndexType *pdest ; /* destination pointer */ 1714 IndexType j ; /* counter */ 1715 IndexType r ; /* a row index */ 1716 IndexType c ; /* a column index */ 1717 IndexType length ; /* length of a row or column */ 1718 1719 /* === Defragment the columns =========================================== */ 1720 1721 pdest = &A[0] ; 1722 for (c = 0 ; c < n_col ; c++) 1723 { 1724 if (COL_IS_ALIVE (c)) 1725 { 1726 psrc = &A [Col [c].start] ; 1727 1728 /* move and compact the column */ 1729 COLAMD_ASSERT (pdest <= psrc) ; 1730 Col [c].start = (IndexType) (pdest - &A [0]) ; 1731 length = Col [c].length ; 1732 for (j = 0 ; j < length ; j++) 1733 { 1734 r = *psrc++ ; 1735 if (ROW_IS_ALIVE (r)) 1736 { 1737 *pdest++ = r ; 1738 } 1739 } 1740 Col [c].length = (IndexType) (pdest - &A [Col [c].start]) ; 1741 } 1742 } 1743 1744 /* === Prepare to defragment the rows =================================== */ 1745 1746 for (r = 0 ; r < n_row ; r++) 1747 { 1748 if (ROW_IS_ALIVE (r)) 1749 { 1750 if (Row [r].length == 0) 1751 { 1752 /* this row is of zero length. cannot compact it, so kill it */ 1753 COLAMD_DEBUG3 (("Defrag row kill\n")) ; 1754 KILL_ROW (r) ; 1755 } 1756 else 1757 { 1758 /* save first column index in Row [r].shared2.first_column */ 1759 psrc = &A [Row [r].start] ; 1760 Row [r].shared2.first_column = *psrc ; 1761 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; 1762 /* flag the start of the row with the one's complement of row */ 1763 *psrc = ONES_COMPLEMENT (r) ; 1764 1765 } 1766 } 1767 } 1768 1769 /* === Defragment the rows ============================================== */ 1770 1771 psrc = pdest ; 1772 while (psrc < pfree) 1773 { 1774 /* find a negative number ... the start of a row */ 1775 if (*psrc++ < 0) 1776 { 1777 psrc-- ; 1778 /* get the row index */ 1779 r = ONES_COMPLEMENT (*psrc) ; 1780 COLAMD_ASSERT (r >= 0 && r < n_row) ; 1781 /* restore first column index */ 1782 *psrc = Row [r].shared2.first_column ; 1783 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; 1784 1785 /* move and compact the row */ 1786 COLAMD_ASSERT (pdest <= psrc) ; 1787 Row [r].start = (IndexType) (pdest - &A [0]) ; 1788 length = Row [r].length ; 1789 for (j = 0 ; j < length ; j++) 1790 { 1791 c = *psrc++ ; 1792 if (COL_IS_ALIVE (c)) 1793 { 1794 *pdest++ = c ; 1795 } 1796 } 1797 Row [r].length = (IndexType) (pdest - &A [Row [r].start]) ; 1798 1799 } 1800 } 1801 /* ensure we found all the rows */ 1802 COLAMD_ASSERT (debug_rows == 0) ; 1803 1804 /* === Return the new value of pfree ==================================== */ 1805 1806 return ((IndexType) (pdest - &A [0])) ; 1807} 1808 1809 1810/* ========================================================================== */ 1811/* === clear_mark =========================================================== */ 1812/* ========================================================================== */ 1813 1814/* 1815 Clears the Row [].shared2.mark array, and returns the new tag_mark. 1816 Return value is the new tag_mark. Not user-callable. 1817*/ 1818template <typename IndexType> 1819static inline IndexType clear_mark /* return the new value for tag_mark */ 1820 ( 1821 /* === Parameters ======================================================= */ 1822 1823 IndexType n_row, /* number of rows in A */ 1824 Colamd_Row<IndexType> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ 1825 ) 1826{ 1827 /* === Local variables ================================================== */ 1828 1829 IndexType r ; 1830 1831 for (r = 0 ; r < n_row ; r++) 1832 { 1833 if (ROW_IS_ALIVE (r)) 1834 { 1835 Row [r].shared2.mark = 0 ; 1836 } 1837 } 1838 return (1) ; 1839} 1840 1841 1842} // namespace internal 1843#endif 1844