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