![]() |
Eigen
3.3.3
|
00001 // // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 // This file is modified from the colamd/symamd library. The copyright is below 00011 00012 // The authors of the code itself are Stefan I. Larimore and Timothy A. 00013 // Davis (davis@cise.ufl.edu), University of Florida. The algorithm was 00014 // developed in collaboration with John Gilbert, Xerox PARC, and Esmond 00015 // Ng, Oak Ridge National Laboratory. 00016 // 00017 // Date: 00018 // 00019 // September 8, 2003. Version 2.3. 00020 // 00021 // Acknowledgements: 00022 // 00023 // This work was supported by the National Science Foundation, under 00024 // grants DMS-9504974 and DMS-9803599. 00025 // 00026 // Notice: 00027 // 00028 // Copyright (c) 1998-2003 by the University of Florida. 00029 // All Rights Reserved. 00030 // 00031 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 00032 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 00033 // 00034 // Permission is hereby granted to use, copy, modify, and/or distribute 00035 // this program, provided that the Copyright, this License, and the 00036 // Availability of the original version is retained on all copies and made 00037 // accessible to the end-user of any code or package that includes COLAMD 00038 // or any modified version of COLAMD. 00039 // 00040 // Availability: 00041 // 00042 // The colamd/symamd library is available at 00043 // 00044 // http://www.suitesparse.com 00045 00046 00047 #ifndef EIGEN_COLAMD_H 00048 #define EIGEN_COLAMD_H 00049 00050 namespace internal { 00051 /* Ensure that debugging is turned off: */ 00052 #ifndef COLAMD_NDEBUG 00053 #define COLAMD_NDEBUG 00054 #endif /* NDEBUG */ 00055 /* ========================================================================== */ 00056 /* === Knob and statistics definitions ====================================== */ 00057 /* ========================================================================== */ 00058 00059 /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ 00060 #define COLAMD_KNOBS 20 00061 00062 /* number of output statistics. Only stats [0..6] are currently used. */ 00063 #define COLAMD_STATS 20 00064 00065 /* knobs [0] and stats [0]: dense row knob and output statistic. */ 00066 #define COLAMD_DENSE_ROW 0 00067 00068 /* knobs [1] and stats [1]: dense column knob and output statistic. */ 00069 #define COLAMD_DENSE_COL 1 00070 00071 /* stats [2]: memory defragmentation count output statistic */ 00072 #define COLAMD_DEFRAG_COUNT 2 00073 00074 /* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ 00075 #define COLAMD_STATUS 3 00076 00077 /* stats [4..6]: error info, or info on jumbled columns */ 00078 #define COLAMD_INFO1 4 00079 #define COLAMD_INFO2 5 00080 #define COLAMD_INFO3 6 00081 00082 /* error codes returned in stats [3]: */ 00083 #define COLAMD_OK (0) 00084 #define COLAMD_OK_BUT_JUMBLED (1) 00085 #define COLAMD_ERROR_A_not_present (-1) 00086 #define COLAMD_ERROR_p_not_present (-2) 00087 #define COLAMD_ERROR_nrow_negative (-3) 00088 #define COLAMD_ERROR_ncol_negative (-4) 00089 #define COLAMD_ERROR_nnz_negative (-5) 00090 #define COLAMD_ERROR_p0_nonzero (-6) 00091 #define COLAMD_ERROR_A_too_small (-7) 00092 #define COLAMD_ERROR_col_length_negative (-8) 00093 #define COLAMD_ERROR_row_index_out_of_bounds (-9) 00094 #define COLAMD_ERROR_out_of_memory (-10) 00095 #define COLAMD_ERROR_internal_error (-999) 00096 00097 /* ========================================================================== */ 00098 /* === Definitions ========================================================== */ 00099 /* ========================================================================== */ 00100 00101 #define ONES_COMPLEMENT(r) (-(r)-1) 00102 00103 /* -------------------------------------------------------------------------- */ 00104 00105 #define COLAMD_EMPTY (-1) 00106 00107 /* Row and column status */ 00108 #define ALIVE (0) 00109 #define DEAD (-1) 00110 00111 /* Column status */ 00112 #define DEAD_PRINCIPAL (-1) 00113 #define DEAD_NON_PRINCIPAL (-2) 00114 00115 /* Macros for row and column status update and checking. */ 00116 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) 00117 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) 00118 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) 00119 #define COL_IS_DEAD(c) (Col [c].start < ALIVE) 00120 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) 00121 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) 00122 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } 00123 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } 00124 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } 00125 00126 /* ========================================================================== */ 00127 /* === Colamd reporting mechanism =========================================== */ 00128 /* ========================================================================== */ 00129 00130 // == Row and Column structures == 00131 template <typename IndexType> 00132 struct colamd_col 00133 { 00134 IndexType start ; /* index for A of first row in this column, or DEAD */ 00135 /* if column is dead */ 00136 IndexType length ; /* number of rows in this column */ 00137 union 00138 { 00139 IndexType thickness ; /* number of original columns represented by this */ 00140 /* col, if the column is alive */ 00141 IndexType parent ; /* parent in parent tree super-column structure, if */ 00142 /* the column is dead */ 00143 } shared1 ; 00144 union 00145 { 00146 IndexType score ; /* the score used to maintain heap, if col is alive */ 00147 IndexType order ; /* pivot ordering of this column, if col is dead */ 00148 } shared2 ; 00149 union 00150 { 00151 IndexType headhash ; /* head of a hash bucket, if col is at the head of */ 00152 /* a degree list */ 00153 IndexType hash ; /* hash value, if col is not in a degree list */ 00154 IndexType prev ; /* previous column in degree list, if col is in a */ 00155 /* degree list (but not at the head of a degree list) */ 00156 } shared3 ; 00157 union 00158 { 00159 IndexType degree_next ; /* next column, if col is in a degree list */ 00160 IndexType hash_next ; /* next column, if col is in a hash list */ 00161 } shared4 ; 00162 00163 }; 00164 00165 template <typename IndexType> 00166 struct Colamd_Row 00167 { 00168 IndexType start ; /* index for A of first col in this row */ 00169 IndexType length ; /* number of principal columns in this row */ 00170 union 00171 { 00172 IndexType degree ; /* number of principal & non-principal columns in row */ 00173 IndexType p ; /* used as a row pointer in init_rows_cols () */ 00174 } shared1 ; 00175 union 00176 { 00177 IndexType mark ; /* for computing set differences and marking dead rows*/ 00178 IndexType first_column ;/* first column in row (used in garbage collection) */ 00179 } shared2 ; 00180 00181 }; 00182 00183 /* ========================================================================== */ 00184 /* === Colamd recommended memory size ======================================= */ 00185 /* ========================================================================== */ 00186 00187 /* 00188 The recommended length Alen of the array A passed to colamd is given by 00189 the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any 00190 argument is negative. 2*nnz space is required for the row and column 00191 indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is 00192 required for the Col and Row arrays, respectively, which are internal to 00193 colamd. An additional n_col space is the minimal amount of "elbow room", 00194 and nnz/5 more space is recommended for run time efficiency. 00195 00196 This macro is not needed when using symamd. 00197 00198 Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid 00199 gcc -pedantic warning messages. 00200 */ 00201 template <typename IndexType> 00202 inline IndexType colamd_c(IndexType n_col) 00203 { return IndexType( ((n_col) + 1) * sizeof (colamd_col<IndexType>) / sizeof (IndexType) ) ; } 00204 00205 template <typename IndexType> 00206 inline IndexType colamd_r(IndexType n_row) 00207 { return IndexType(((n_row) + 1) * sizeof (Colamd_Row<IndexType>) / sizeof (IndexType)); } 00208 00209 // Prototypes of non-user callable routines 00210 template <typename IndexType> 00211 static 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] ); 00212 00213 template <typename IndexType> 00214 static 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); 00215 00216 template <typename IndexType> 00217 static 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); 00218 00219 template <typename IndexType> 00220 static void order_children (IndexType n_col, colamd_col<IndexType> Col [], IndexType p []); 00221 00222 template <typename IndexType> 00223 static void detect_super_cols (colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType row_start, IndexType row_length ) ; 00224 00225 template <typename IndexType> 00226 static IndexType garbage_collection (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType *pfree) ; 00227 00228 template <typename IndexType> 00229 static inline IndexType clear_mark (IndexType n_row, Colamd_Row<IndexType> Row [] ) ; 00230 00231 /* === No debugging ========================================================= */ 00232 00233 #define COLAMD_DEBUG0(params) ; 00234 #define COLAMD_DEBUG1(params) ; 00235 #define COLAMD_DEBUG2(params) ; 00236 #define COLAMD_DEBUG3(params) ; 00237 #define COLAMD_DEBUG4(params) ; 00238 00239 #define COLAMD_ASSERT(expression) ((void) 0) 00240 00241 00256 template <typename IndexType> 00257 inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col) 00258 { 00259 if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) 00260 return (-1); 00261 else 00262 return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5)); 00263 } 00264 00286 static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS]) 00287 { 00288 /* === Local variables ================================================== */ 00289 00290 int i ; 00291 00292 if (!knobs) 00293 { 00294 return ; /* no knobs to initialize */ 00295 } 00296 for (i = 0 ; i < COLAMD_KNOBS ; i++) 00297 { 00298 knobs [i] = 0 ; 00299 } 00300 knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ 00301 knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ 00302 } 00303 00321 template <typename IndexType> 00322 static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS]) 00323 { 00324 /* === Local variables ================================================== */ 00325 00326 IndexType i ; /* loop index */ 00327 IndexType nnz ; /* nonzeros in A */ 00328 IndexType Row_size ; /* size of Row [], in integers */ 00329 IndexType Col_size ; /* size of Col [], in integers */ 00330 IndexType need ; /* minimum required length of A */ 00331 Colamd_Row<IndexType> *Row ; /* pointer into A of Row [0..n_row] array */ 00332 colamd_col<IndexType> *Col ; /* pointer into A of Col [0..n_col] array */ 00333 IndexType n_col2 ; /* number of non-dense, non-empty columns */ 00334 IndexType n_row2 ; /* number of non-dense, non-empty rows */ 00335 IndexType ngarbage ; /* number of garbage collections performed */ 00336 IndexType max_deg ; /* maximum row degree */ 00337 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ 00338 00339 00340 /* === Check the input arguments ======================================== */ 00341 00342 if (!stats) 00343 { 00344 COLAMD_DEBUG0 (("colamd: stats not present\n")) ; 00345 return (false) ; 00346 } 00347 for (i = 0 ; i < COLAMD_STATS ; i++) 00348 { 00349 stats [i] = 0 ; 00350 } 00351 stats [COLAMD_STATUS] = COLAMD_OK ; 00352 stats [COLAMD_INFO1] = -1 ; 00353 stats [COLAMD_INFO2] = -1 ; 00354 00355 if (!A) /* A is not present */ 00356 { 00357 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; 00358 COLAMD_DEBUG0 (("colamd: A not present\n")) ; 00359 return (false) ; 00360 } 00361 00362 if (!p) /* p is not present */ 00363 { 00364 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; 00365 COLAMD_DEBUG0 (("colamd: p not present\n")) ; 00366 return (false) ; 00367 } 00368 00369 if (n_row < 0) /* n_row must be >= 0 */ 00370 { 00371 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; 00372 stats [COLAMD_INFO1] = n_row ; 00373 COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; 00374 return (false) ; 00375 } 00376 00377 if (n_col < 0) /* n_col must be >= 0 */ 00378 { 00379 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; 00380 stats [COLAMD_INFO1] = n_col ; 00381 COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; 00382 return (false) ; 00383 } 00384 00385 nnz = p [n_col] ; 00386 if (nnz < 0) /* nnz must be >= 0 */ 00387 { 00388 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; 00389 stats [COLAMD_INFO1] = nnz ; 00390 COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; 00391 return (false) ; 00392 } 00393 00394 if (p [0] != 0) 00395 { 00396 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; 00397 stats [COLAMD_INFO1] = p [0] ; 00398 COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; 00399 return (false) ; 00400 } 00401 00402 /* === If no knobs, set default knobs =================================== */ 00403 00404 if (!knobs) 00405 { 00406 colamd_set_defaults (default_knobs) ; 00407 knobs = default_knobs ; 00408 } 00409 00410 /* === Allocate the Row and Col arrays from array A ===================== */ 00411 00412 Col_size = colamd_c (n_col) ; 00413 Row_size = colamd_r (n_row) ; 00414 need = 2*nnz + n_col + Col_size + Row_size ; 00415 00416 if (need > Alen) 00417 { 00418 /* not enough space in array A to perform the ordering */ 00419 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; 00420 stats [COLAMD_INFO1] = need ; 00421 stats [COLAMD_INFO2] = Alen ; 00422 COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); 00423 return (false) ; 00424 } 00425 00426 Alen -= Col_size + Row_size ; 00427 Col = (colamd_col<IndexType> *) &A [Alen] ; 00428 Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ; 00429 00430 /* === Construct the row and column data structures ===================== */ 00431 00432 if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) 00433 { 00434 /* input matrix is invalid */ 00435 COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ; 00436 return (false) ; 00437 } 00438 00439 /* === Initialize scores, kill dense rows/columns ======================= */ 00440 00441 Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs, 00442 &n_row2, &n_col2, &max_deg) ; 00443 00444 /* === Order the supercolumns =========================================== */ 00445 00446 ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p, 00447 n_col2, max_deg, 2*nnz) ; 00448 00449 /* === Order the non-principal columns ================================== */ 00450 00451 Eigen::internal::order_children (n_col, Col, p) ; 00452 00453 /* === Return statistics in stats ======================================= */ 00454 00455 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; 00456 stats [COLAMD_DENSE_COL] = n_col - n_col2 ; 00457 stats [COLAMD_DEFRAG_COUNT] = ngarbage ; 00458 COLAMD_DEBUG0 (("colamd: done.\n")) ; 00459 return (true) ; 00460 } 00461 00462 /* ========================================================================== */ 00463 /* === NON-USER-CALLABLE ROUTINES: ========================================== */ 00464 /* ========================================================================== */ 00465 00466 /* There are no user-callable routines beyond this point in the file */ 00467 00468 00469 /* ========================================================================== */ 00470 /* === init_rows_cols ======================================================= */ 00471 /* ========================================================================== */ 00472 00473 /* 00474 Takes the column form of the matrix in A and creates the row form of the 00475 matrix. Also, row and column attributes are stored in the Col and Row 00476 structs. If the columns are un-sorted or contain duplicate row indices, 00477 this routine will also sort and remove duplicate row indices from the 00478 column form of the matrix. Returns false if the matrix is invalid, 00479 true otherwise. Not user-callable. 00480 */ 00481 template <typename IndexType> 00482 static IndexType init_rows_cols /* returns true if OK, or false otherwise */ 00483 ( 00484 /* === Parameters ======================================================= */ 00485 00486 IndexType n_row, /* number of rows of A */ 00487 IndexType n_col, /* number of columns of A */ 00488 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 00489 colamd_col<IndexType> Col [], /* of size n_col+1 */ 00490 IndexType A [], /* row indices of A, of size Alen */ 00491 IndexType p [], /* pointers to columns in A, of size n_col+1 */ 00492 IndexType stats [COLAMD_STATS] /* colamd statistics */ 00493 ) 00494 { 00495 /* === Local variables ================================================== */ 00496 00497 IndexType col ; /* a column index */ 00498 IndexType row ; /* a row index */ 00499 IndexType *cp ; /* a column pointer */ 00500 IndexType *cp_end ; /* a pointer to the end of a column */ 00501 IndexType *rp ; /* a row pointer */ 00502 IndexType *rp_end ; /* a pointer to the end of a row */ 00503 IndexType last_row ; /* previous row */ 00504 00505 /* === Initialize columns, and check column pointers ==================== */ 00506 00507 for (col = 0 ; col < n_col ; col++) 00508 { 00509 Col [col].start = p [col] ; 00510 Col [col].length = p [col+1] - p [col] ; 00511 00512 if ((Col [col].length) < 0) // extra parentheses to work-around gcc bug 10200 00513 { 00514 /* column pointers must be non-decreasing */ 00515 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; 00516 stats [COLAMD_INFO1] = col ; 00517 stats [COLAMD_INFO2] = Col [col].length ; 00518 COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; 00519 return (false) ; 00520 } 00521 00522 Col [col].shared1.thickness = 1 ; 00523 Col [col].shared2.score = 0 ; 00524 Col [col].shared3.prev = COLAMD_EMPTY ; 00525 Col [col].shared4.degree_next = COLAMD_EMPTY ; 00526 } 00527 00528 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ 00529 00530 /* === Scan columns, compute row degrees, and check row indices ========= */ 00531 00532 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ 00533 00534 for (row = 0 ; row < n_row ; row++) 00535 { 00536 Row [row].length = 0 ; 00537 Row [row].shared2.mark = -1 ; 00538 } 00539 00540 for (col = 0 ; col < n_col ; col++) 00541 { 00542 last_row = -1 ; 00543 00544 cp = &A [p [col]] ; 00545 cp_end = &A [p [col+1]] ; 00546 00547 while (cp < cp_end) 00548 { 00549 row = *cp++ ; 00550 00551 /* make sure row indices within range */ 00552 if (row < 0 || row >= n_row) 00553 { 00554 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; 00555 stats [COLAMD_INFO1] = col ; 00556 stats [COLAMD_INFO2] = row ; 00557 stats [COLAMD_INFO3] = n_row ; 00558 COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; 00559 return (false) ; 00560 } 00561 00562 if (row <= last_row || Row [row].shared2.mark == col) 00563 { 00564 /* row index are unsorted or repeated (or both), thus col */ 00565 /* is jumbled. This is a notice, not an error condition. */ 00566 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; 00567 stats [COLAMD_INFO1] = col ; 00568 stats [COLAMD_INFO2] = row ; 00569 (stats [COLAMD_INFO3]) ++ ; 00570 COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); 00571 } 00572 00573 if (Row [row].shared2.mark != col) 00574 { 00575 Row [row].length++ ; 00576 } 00577 else 00578 { 00579 /* this is a repeated entry in the column, */ 00580 /* it will be removed */ 00581 Col [col].length-- ; 00582 } 00583 00584 /* mark the row as having been seen in this column */ 00585 Row [row].shared2.mark = col ; 00586 00587 last_row = row ; 00588 } 00589 } 00590 00591 /* === Compute row pointers ============================================= */ 00592 00593 /* row form of the matrix starts directly after the column */ 00594 /* form of matrix in A */ 00595 Row [0].start = p [n_col] ; 00596 Row [0].shared1.p = Row [0].start ; 00597 Row [0].shared2.mark = -1 ; 00598 for (row = 1 ; row < n_row ; row++) 00599 { 00600 Row [row].start = Row [row-1].start + Row [row-1].length ; 00601 Row [row].shared1.p = Row [row].start ; 00602 Row [row].shared2.mark = -1 ; 00603 } 00604 00605 /* === Create row form ================================================== */ 00606 00607 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 00608 { 00609 /* if cols jumbled, watch for repeated row indices */ 00610 for (col = 0 ; col < n_col ; col++) 00611 { 00612 cp = &A [p [col]] ; 00613 cp_end = &A [p [col+1]] ; 00614 while (cp < cp_end) 00615 { 00616 row = *cp++ ; 00617 if (Row [row].shared2.mark != col) 00618 { 00619 A [(Row [row].shared1.p)++] = col ; 00620 Row [row].shared2.mark = col ; 00621 } 00622 } 00623 } 00624 } 00625 else 00626 { 00627 /* if cols not jumbled, we don't need the mark (this is faster) */ 00628 for (col = 0 ; col < n_col ; col++) 00629 { 00630 cp = &A [p [col]] ; 00631 cp_end = &A [p [col+1]] ; 00632 while (cp < cp_end) 00633 { 00634 A [(Row [*cp++].shared1.p)++] = col ; 00635 } 00636 } 00637 } 00638 00639 /* === Clear the row marks and set row degrees ========================== */ 00640 00641 for (row = 0 ; row < n_row ; row++) 00642 { 00643 Row [row].shared2.mark = 0 ; 00644 Row [row].shared1.degree = Row [row].length ; 00645 } 00646 00647 /* === See if we need to re-create columns ============================== */ 00648 00649 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) 00650 { 00651 COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; 00652 00653 00654 /* === Compute col pointers ========================================= */ 00655 00656 /* col form of the matrix starts at A [0]. */ 00657 /* Note, we may have a gap between the col form and the row */ 00658 /* form if there were duplicate entries, if so, it will be */ 00659 /* removed upon the first garbage collection */ 00660 Col [0].start = 0 ; 00661 p [0] = Col [0].start ; 00662 for (col = 1 ; col < n_col ; col++) 00663 { 00664 /* note that the lengths here are for pruned columns, i.e. */ 00665 /* no duplicate row indices will exist for these columns */ 00666 Col [col].start = Col [col-1].start + Col [col-1].length ; 00667 p [col] = Col [col].start ; 00668 } 00669 00670 /* === Re-create col form =========================================== */ 00671 00672 for (row = 0 ; row < n_row ; row++) 00673 { 00674 rp = &A [Row [row].start] ; 00675 rp_end = rp + Row [row].length ; 00676 while (rp < rp_end) 00677 { 00678 A [(p [*rp++])++] = row ; 00679 } 00680 } 00681 } 00682 00683 /* === Done. Matrix is not (or no longer) jumbled ====================== */ 00684 00685 return (true) ; 00686 } 00687 00688 00689 /* ========================================================================== */ 00690 /* === init_scoring ========================================================= */ 00691 /* ========================================================================== */ 00692 00693 /* 00694 Kills dense or empty columns and rows, calculates an initial score for 00695 each column, and places all columns in the degree lists. Not user-callable. 00696 */ 00697 template <typename IndexType> 00698 static void init_scoring 00699 ( 00700 /* === Parameters ======================================================= */ 00701 00702 IndexType n_row, /* number of rows of A */ 00703 IndexType n_col, /* number of columns of A */ 00704 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 00705 colamd_col<IndexType> Col [], /* of size n_col+1 */ 00706 IndexType A [], /* column form and row form of A */ 00707 IndexType head [], /* of size n_col+1 */ 00708 double knobs [COLAMD_KNOBS],/* parameters */ 00709 IndexType *p_n_row2, /* number of non-dense, non-empty rows */ 00710 IndexType *p_n_col2, /* number of non-dense, non-empty columns */ 00711 IndexType *p_max_deg /* maximum row degree */ 00712 ) 00713 { 00714 /* === Local variables ================================================== */ 00715 00716 IndexType c ; /* a column index */ 00717 IndexType r, row ; /* a row index */ 00718 IndexType *cp ; /* a column pointer */ 00719 IndexType deg ; /* degree of a row or column */ 00720 IndexType *cp_end ; /* a pointer to the end of a column */ 00721 IndexType *new_cp ; /* new column pointer */ 00722 IndexType col_length ; /* length of pruned column */ 00723 IndexType score ; /* current column score */ 00724 IndexType n_col2 ; /* number of non-dense, non-empty columns */ 00725 IndexType n_row2 ; /* number of non-dense, non-empty rows */ 00726 IndexType dense_row_count ; /* remove rows with more entries than this */ 00727 IndexType dense_col_count ; /* remove cols with more entries than this */ 00728 IndexType min_score ; /* smallest column score */ 00729 IndexType max_deg ; /* maximum row degree */ 00730 IndexType next_col ; /* Used to add to degree list.*/ 00731 00732 00733 /* === Extract knobs ==================================================== */ 00734 00735 dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ; 00736 dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ; 00737 COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; 00738 max_deg = 0 ; 00739 n_col2 = n_col ; 00740 n_row2 = n_row ; 00741 00742 /* === Kill empty columns =============================================== */ 00743 00744 /* Put the empty columns at the end in their natural order, so that LU */ 00745 /* factorization can proceed as far as possible. */ 00746 for (c = n_col-1 ; c >= 0 ; c--) 00747 { 00748 deg = Col [c].length ; 00749 if (deg == 0) 00750 { 00751 /* this is a empty column, kill and order it last */ 00752 Col [c].shared2.order = --n_col2 ; 00753 KILL_PRINCIPAL_COL (c) ; 00754 } 00755 } 00756 COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; 00757 00758 /* === Kill dense columns =============================================== */ 00759 00760 /* Put the dense columns at the end, in their natural order */ 00761 for (c = n_col-1 ; c >= 0 ; c--) 00762 { 00763 /* skip any dead columns */ 00764 if (COL_IS_DEAD (c)) 00765 { 00766 continue ; 00767 } 00768 deg = Col [c].length ; 00769 if (deg > dense_col_count) 00770 { 00771 /* this is a dense column, kill and order it last */ 00772 Col [c].shared2.order = --n_col2 ; 00773 /* decrement the row degrees */ 00774 cp = &A [Col [c].start] ; 00775 cp_end = cp + Col [c].length ; 00776 while (cp < cp_end) 00777 { 00778 Row [*cp++].shared1.degree-- ; 00779 } 00780 KILL_PRINCIPAL_COL (c) ; 00781 } 00782 } 00783 COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; 00784 00785 /* === Kill dense and empty rows ======================================== */ 00786 00787 for (r = 0 ; r < n_row ; r++) 00788 { 00789 deg = Row [r].shared1.degree ; 00790 COLAMD_ASSERT (deg >= 0 && deg <= n_col) ; 00791 if (deg > dense_row_count || deg == 0) 00792 { 00793 /* kill a dense or empty row */ 00794 KILL_ROW (r) ; 00795 --n_row2 ; 00796 } 00797 else 00798 { 00799 /* keep track of max degree of remaining rows */ 00800 max_deg = numext::maxi(max_deg, deg) ; 00801 } 00802 } 00803 COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; 00804 00805 /* === Compute initial column scores ==================================== */ 00806 00807 /* At this point the row degrees are accurate. They reflect the number */ 00808 /* of "live" (non-dense) columns in each row. No empty rows exist. */ 00809 /* Some "live" columns may contain only dead rows, however. These are */ 00810 /* pruned in the code below. */ 00811 00812 /* now find the initial matlab score for each column */ 00813 for (c = n_col-1 ; c >= 0 ; c--) 00814 { 00815 /* skip dead column */ 00816 if (COL_IS_DEAD (c)) 00817 { 00818 continue ; 00819 } 00820 score = 0 ; 00821 cp = &A [Col [c].start] ; 00822 new_cp = cp ; 00823 cp_end = cp + Col [c].length ; 00824 while (cp < cp_end) 00825 { 00826 /* get a row */ 00827 row = *cp++ ; 00828 /* skip if dead */ 00829 if (ROW_IS_DEAD (row)) 00830 { 00831 continue ; 00832 } 00833 /* compact the column */ 00834 *new_cp++ = row ; 00835 /* add row's external degree */ 00836 score += Row [row].shared1.degree - 1 ; 00837 /* guard against integer overflow */ 00838 score = numext::mini(score, n_col) ; 00839 } 00840 /* determine pruned column length */ 00841 col_length = (IndexType) (new_cp - &A [Col [c].start]) ; 00842 if (col_length == 0) 00843 { 00844 /* a newly-made null column (all rows in this col are "dense" */ 00845 /* and have already been killed) */ 00846 COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ; 00847 Col [c].shared2.order = --n_col2 ; 00848 KILL_PRINCIPAL_COL (c) ; 00849 } 00850 else 00851 { 00852 /* set column length and set score */ 00853 COLAMD_ASSERT (score >= 0) ; 00854 COLAMD_ASSERT (score <= n_col) ; 00855 Col [c].length = col_length ; 00856 Col [c].shared2.score = score ; 00857 } 00858 } 00859 COLAMD_DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", 00860 n_col-n_col2)) ; 00861 00862 /* At this point, all empty rows and columns are dead. All live columns */ 00863 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ 00864 /* yet). Rows may contain dead columns, but all live rows contain at */ 00865 /* least one live column. */ 00866 00867 /* === Initialize degree lists ========================================== */ 00868 00869 00870 /* clear the hash buckets */ 00871 for (c = 0 ; c <= n_col ; c++) 00872 { 00873 head [c] = COLAMD_EMPTY ; 00874 } 00875 min_score = n_col ; 00876 /* place in reverse order, so low column indices are at the front */ 00877 /* of the lists. This is to encourage natural tie-breaking */ 00878 for (c = n_col-1 ; c >= 0 ; c--) 00879 { 00880 /* only add principal columns to degree lists */ 00881 if (COL_IS_ALIVE (c)) 00882 { 00883 COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n", 00884 c, Col [c].shared2.score, min_score, n_col)) ; 00885 00886 /* === Add columns score to DList =============================== */ 00887 00888 score = Col [c].shared2.score ; 00889 00890 COLAMD_ASSERT (min_score >= 0) ; 00891 COLAMD_ASSERT (min_score <= n_col) ; 00892 COLAMD_ASSERT (score >= 0) ; 00893 COLAMD_ASSERT (score <= n_col) ; 00894 COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ; 00895 00896 /* now add this column to dList at proper score location */ 00897 next_col = head [score] ; 00898 Col [c].shared3.prev = COLAMD_EMPTY ; 00899 Col [c].shared4.degree_next = next_col ; 00900 00901 /* if there already was a column with the same score, set its */ 00902 /* previous pointer to this new column */ 00903 if (next_col != COLAMD_EMPTY) 00904 { 00905 Col [next_col].shared3.prev = c ; 00906 } 00907 head [score] = c ; 00908 00909 /* see if this score is less than current min */ 00910 min_score = numext::mini(min_score, score) ; 00911 00912 00913 } 00914 } 00915 00916 00917 /* === Return number of remaining columns, and max row degree =========== */ 00918 00919 *p_n_col2 = n_col2 ; 00920 *p_n_row2 = n_row2 ; 00921 *p_max_deg = max_deg ; 00922 } 00923 00924 00925 /* ========================================================================== */ 00926 /* === find_ordering ======================================================== */ 00927 /* ========================================================================== */ 00928 00929 /* 00930 Order the principal columns of the supercolumn form of the matrix 00931 (no supercolumns on input). Uses a minimum approximate column minimum 00932 degree ordering method. Not user-callable. 00933 */ 00934 template <typename IndexType> 00935 static IndexType find_ordering /* return the number of garbage collections */ 00936 ( 00937 /* === Parameters ======================================================= */ 00938 00939 IndexType n_row, /* number of rows of A */ 00940 IndexType n_col, /* number of columns of A */ 00941 IndexType Alen, /* size of A, 2*nnz + n_col or larger */ 00942 Colamd_Row<IndexType> Row [], /* of size n_row+1 */ 00943 colamd_col<IndexType> Col [], /* of size n_col+1 */ 00944 IndexType A [], /* column form and row form of A */ 00945 IndexType head [], /* of size n_col+1 */ 00946 IndexType n_col2, /* Remaining columns to order */ 00947 IndexType max_deg, /* Maximum row degree */ 00948 IndexType pfree /* index of first free slot (2*nnz on entry) */ 00949 ) 00950 { 00951 /* === Local variables ================================================== */ 00952 00953 IndexType k ; /* current pivot ordering step */ 00954 IndexType pivot_col ; /* current pivot column */ 00955 IndexType *cp ; /* a column pointer */ 00956 IndexType *rp ; /* a row pointer */ 00957 IndexType pivot_row ; /* current pivot row */ 00958 IndexType *new_cp ; /* modified column pointer */ 00959 IndexType *new_rp ; /* modified row pointer */ 00960 IndexType pivot_row_start ; /* pointer to start of pivot row */ 00961 IndexType pivot_row_degree ; /* number of columns in pivot row */ 00962 IndexType pivot_row_length ; /* number of supercolumns in pivot row */ 00963 IndexType pivot_col_score ; /* score of pivot column */ 00964 IndexType needed_memory ; /* free space needed for pivot row */ 00965 IndexType *cp_end ; /* pointer to the end of a column */ 00966 IndexType *rp_end ; /* pointer to the end of a row */ 00967 IndexType row ; /* a row index */ 00968 IndexType col ; /* a column index */ 00969 IndexType max_score ; /* maximum possible score */ 00970 IndexType cur_score ; /* score of current column */ 00971 unsigned int hash ; /* hash value for supernode detection */ 00972 IndexType head_column ; /* head of hash bucket */ 00973 IndexType first_col ; /* first column in hash bucket */ 00974 IndexType tag_mark ; /* marker value for mark array */ 00975 IndexType row_mark ; /* Row [row].shared2.mark */ 00976 IndexType set_difference ; /* set difference size of row with pivot row */ 00977 IndexType min_score ; /* smallest column score */ 00978 IndexType col_thickness ; /* "thickness" (no. of columns in a supercol) */ 00979 IndexType max_mark ; /* maximum value of tag_mark */ 00980 IndexType pivot_col_thickness ; /* number of columns represented by pivot col */ 00981 IndexType prev_col ; /* Used by Dlist operations. */ 00982 IndexType next_col ; /* Used by Dlist operations. */ 00983 IndexType ngarbage ; /* number of garbage collections performed */ 00984 00985 00986 /* === Initialization and clear mark ==================================== */ 00987 00988 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */ 00989 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 00990 min_score = 0 ; 00991 ngarbage = 0 ; 00992 COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; 00993 00994 /* === Order the columns ================================================ */ 00995 00996 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) 00997 { 00998 00999 /* === Select pivot column, and order it ============================ */ 01000 01001 /* make sure degree list isn't empty */ 01002 COLAMD_ASSERT (min_score >= 0) ; 01003 COLAMD_ASSERT (min_score <= n_col) ; 01004 COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ; 01005 01006 /* get pivot column from head of minimum degree list */ 01007 while (head [min_score] == COLAMD_EMPTY && min_score < n_col) 01008 { 01009 min_score++ ; 01010 } 01011 pivot_col = head [min_score] ; 01012 COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; 01013 next_col = Col [pivot_col].shared4.degree_next ; 01014 head [min_score] = next_col ; 01015 if (next_col != COLAMD_EMPTY) 01016 { 01017 Col [next_col].shared3.prev = COLAMD_EMPTY ; 01018 } 01019 01020 COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ; 01021 COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ; 01022 01023 /* remember score for defrag check */ 01024 pivot_col_score = Col [pivot_col].shared2.score ; 01025 01026 /* the pivot column is the kth column in the pivot order */ 01027 Col [pivot_col].shared2.order = k ; 01028 01029 /* increment order count by column thickness */ 01030 pivot_col_thickness = Col [pivot_col].shared1.thickness ; 01031 k += pivot_col_thickness ; 01032 COLAMD_ASSERT (pivot_col_thickness > 0) ; 01033 01034 /* === Garbage_collection, if necessary ============================= */ 01035 01036 needed_memory = numext::mini(pivot_col_score, n_col - k) ; 01037 if (pfree + needed_memory >= Alen) 01038 { 01039 pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; 01040 ngarbage++ ; 01041 /* after garbage collection we will have enough */ 01042 COLAMD_ASSERT (pfree + needed_memory < Alen) ; 01043 /* garbage collection has wiped out the Row[].shared2.mark array */ 01044 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 01045 01046 } 01047 01048 /* === Compute pivot row pattern ==================================== */ 01049 01050 /* get starting location for this new merged row */ 01051 pivot_row_start = pfree ; 01052 01053 /* initialize new row counts to zero */ 01054 pivot_row_degree = 0 ; 01055 01056 /* tag pivot column as having been visited so it isn't included */ 01057 /* in merged pivot row */ 01058 Col [pivot_col].shared1.thickness = -pivot_col_thickness ; 01059 01060 /* pivot row is the union of all rows in the pivot column pattern */ 01061 cp = &A [Col [pivot_col].start] ; 01062 cp_end = cp + Col [pivot_col].length ; 01063 while (cp < cp_end) 01064 { 01065 /* get a row */ 01066 row = *cp++ ; 01067 COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; 01068 /* skip if row is dead */ 01069 if (ROW_IS_DEAD (row)) 01070 { 01071 continue ; 01072 } 01073 rp = &A [Row [row].start] ; 01074 rp_end = rp + Row [row].length ; 01075 while (rp < rp_end) 01076 { 01077 /* get a column */ 01078 col = *rp++ ; 01079 /* add the column, if alive and untagged */ 01080 col_thickness = Col [col].shared1.thickness ; 01081 if (col_thickness > 0 && COL_IS_ALIVE (col)) 01082 { 01083 /* tag column in pivot row */ 01084 Col [col].shared1.thickness = -col_thickness ; 01085 COLAMD_ASSERT (pfree < Alen) ; 01086 /* place column in pivot row */ 01087 A [pfree++] = col ; 01088 pivot_row_degree += col_thickness ; 01089 } 01090 } 01091 } 01092 01093 /* clear tag on pivot column */ 01094 Col [pivot_col].shared1.thickness = pivot_col_thickness ; 01095 max_deg = numext::maxi(max_deg, pivot_row_degree) ; 01096 01097 01098 /* === Kill all rows used to construct pivot row ==================== */ 01099 01100 /* also kill pivot row, temporarily */ 01101 cp = &A [Col [pivot_col].start] ; 01102 cp_end = cp + Col [pivot_col].length ; 01103 while (cp < cp_end) 01104 { 01105 /* may be killing an already dead row */ 01106 row = *cp++ ; 01107 COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ; 01108 KILL_ROW (row) ; 01109 } 01110 01111 /* === Select a row index to use as the new pivot row =============== */ 01112 01113 pivot_row_length = pfree - pivot_row_start ; 01114 if (pivot_row_length > 0) 01115 { 01116 /* pick the "pivot" row arbitrarily (first row in col) */ 01117 pivot_row = A [Col [pivot_col].start] ; 01118 COLAMD_DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; 01119 } 01120 else 01121 { 01122 /* there is no pivot row, since it is of zero length */ 01123 pivot_row = COLAMD_EMPTY ; 01124 COLAMD_ASSERT (pivot_row_length == 0) ; 01125 } 01126 COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; 01127 01128 /* === Approximate degree computation =============================== */ 01129 01130 /* Here begins the computation of the approximate degree. The column */ 01131 /* score is the sum of the pivot row "length", plus the size of the */ 01132 /* set differences of each row in the column minus the pattern of the */ 01133 /* pivot row itself. The column ("thickness") itself is also */ 01134 /* excluded from the column score (we thus use an approximate */ 01135 /* external degree). */ 01136 01137 /* The time taken by the following code (compute set differences, and */ 01138 /* add them up) is proportional to the size of the data structure */ 01139 /* being scanned - that is, the sum of the sizes of each column in */ 01140 /* the pivot row. Thus, the amortized time to compute a column score */ 01141 /* is proportional to the size of that column (where size, in this */ 01142 /* context, is the column "length", or the number of row indices */ 01143 /* in that column). The number of row indices in a column is */ 01144 /* monotonically non-decreasing, from the length of the original */ 01145 /* column on input to colamd. */ 01146 01147 /* === Compute set differences ====================================== */ 01148 01149 COLAMD_DEBUG3 (("** Computing set differences phase. **\n")) ; 01150 01151 /* pivot row is currently dead - it will be revived later. */ 01152 01153 COLAMD_DEBUG3 (("Pivot row: ")) ; 01154 /* for each column in pivot row */ 01155 rp = &A [pivot_row_start] ; 01156 rp_end = rp + pivot_row_length ; 01157 while (rp < rp_end) 01158 { 01159 col = *rp++ ; 01160 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 01161 COLAMD_DEBUG3 (("Col: %d\n", col)) ; 01162 01163 /* clear tags used to construct pivot row pattern */ 01164 col_thickness = -Col [col].shared1.thickness ; 01165 COLAMD_ASSERT (col_thickness > 0) ; 01166 Col [col].shared1.thickness = col_thickness ; 01167 01168 /* === Remove column from degree list =========================== */ 01169 01170 cur_score = Col [col].shared2.score ; 01171 prev_col = Col [col].shared3.prev ; 01172 next_col = Col [col].shared4.degree_next ; 01173 COLAMD_ASSERT (cur_score >= 0) ; 01174 COLAMD_ASSERT (cur_score <= n_col) ; 01175 COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ; 01176 if (prev_col == COLAMD_EMPTY) 01177 { 01178 head [cur_score] = next_col ; 01179 } 01180 else 01181 { 01182 Col [prev_col].shared4.degree_next = next_col ; 01183 } 01184 if (next_col != COLAMD_EMPTY) 01185 { 01186 Col [next_col].shared3.prev = prev_col ; 01187 } 01188 01189 /* === Scan the column ========================================== */ 01190 01191 cp = &A [Col [col].start] ; 01192 cp_end = cp + Col [col].length ; 01193 while (cp < cp_end) 01194 { 01195 /* get a row */ 01196 row = *cp++ ; 01197 row_mark = Row [row].shared2.mark ; 01198 /* skip if dead */ 01199 if (ROW_IS_MARKED_DEAD (row_mark)) 01200 { 01201 continue ; 01202 } 01203 COLAMD_ASSERT (row != pivot_row) ; 01204 set_difference = row_mark - tag_mark ; 01205 /* check if the row has been seen yet */ 01206 if (set_difference < 0) 01207 { 01208 COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ; 01209 set_difference = Row [row].shared1.degree ; 01210 } 01211 /* subtract column thickness from this row's set difference */ 01212 set_difference -= col_thickness ; 01213 COLAMD_ASSERT (set_difference >= 0) ; 01214 /* absorb this row if the set difference becomes zero */ 01215 if (set_difference == 0) 01216 { 01217 COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; 01218 KILL_ROW (row) ; 01219 } 01220 else 01221 { 01222 /* save the new mark */ 01223 Row [row].shared2.mark = set_difference + tag_mark ; 01224 } 01225 } 01226 } 01227 01228 01229 /* === Add up set differences for each column ======================= */ 01230 01231 COLAMD_DEBUG3 (("** Adding set differences phase. **\n")) ; 01232 01233 /* for each column in pivot row */ 01234 rp = &A [pivot_row_start] ; 01235 rp_end = rp + pivot_row_length ; 01236 while (rp < rp_end) 01237 { 01238 /* get a column */ 01239 col = *rp++ ; 01240 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; 01241 hash = 0 ; 01242 cur_score = 0 ; 01243 cp = &A [Col [col].start] ; 01244 /* compact the column */ 01245 new_cp = cp ; 01246 cp_end = cp + Col [col].length ; 01247 01248 COLAMD_DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; 01249 01250 while (cp < cp_end) 01251 { 01252 /* get a row */ 01253 row = *cp++ ; 01254 COLAMD_ASSERT(row >= 0 && row < n_row) ; 01255 row_mark = Row [row].shared2.mark ; 01256 /* skip if dead */ 01257 if (ROW_IS_MARKED_DEAD (row_mark)) 01258 { 01259 continue ; 01260 } 01261 COLAMD_ASSERT (row_mark > tag_mark) ; 01262 /* compact the column */ 01263 *new_cp++ = row ; 01264 /* compute hash function */ 01265 hash += row ; 01266 /* add set difference */ 01267 cur_score += row_mark - tag_mark ; 01268 /* integer overflow... */ 01269 cur_score = numext::mini(cur_score, n_col) ; 01270 } 01271 01272 /* recompute the column's length */ 01273 Col [col].length = (IndexType) (new_cp - &A [Col [col].start]) ; 01274 01275 /* === Further mass elimination ================================= */ 01276 01277 if (Col [col].length == 0) 01278 { 01279 COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ; 01280 /* nothing left but the pivot row in this column */ 01281 KILL_PRINCIPAL_COL (col) ; 01282 pivot_row_degree -= Col [col].shared1.thickness ; 01283 COLAMD_ASSERT (pivot_row_degree >= 0) ; 01284 /* order it */ 01285 Col [col].shared2.order = k ; 01286 /* increment order count by column thickness */ 01287 k += Col [col].shared1.thickness ; 01288 } 01289 else 01290 { 01291 /* === Prepare for supercolumn detection ==================== */ 01292 01293 COLAMD_DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; 01294 01295 /* save score so far */ 01296 Col [col].shared2.score = cur_score ; 01297 01298 /* add column to hash table, for supercolumn detection */ 01299 hash %= n_col + 1 ; 01300 01301 COLAMD_DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; 01302 COLAMD_ASSERT (hash <= n_col) ; 01303 01304 head_column = head [hash] ; 01305 if (head_column > COLAMD_EMPTY) 01306 { 01307 /* degree list "hash" is non-empty, use prev (shared3) of */ 01308 /* first column in degree list as head of hash bucket */ 01309 first_col = Col [head_column].shared3.headhash ; 01310 Col [head_column].shared3.headhash = col ; 01311 } 01312 else 01313 { 01314 /* degree list "hash" is empty, use head as hash bucket */ 01315 first_col = - (head_column + 2) ; 01316 head [hash] = - (col + 2) ; 01317 } 01318 Col [col].shared4.hash_next = first_col ; 01319 01320 /* save hash function in Col [col].shared3.hash */ 01321 Col [col].shared3.hash = (IndexType) hash ; 01322 COLAMD_ASSERT (COL_IS_ALIVE (col)) ; 01323 } 01324 } 01325 01326 /* The approximate external column degree is now computed. */ 01327 01328 /* === Supercolumn detection ======================================== */ 01329 01330 COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ; 01331 01332 Eigen::internal::detect_super_cols (Col, A, head, pivot_row_start, pivot_row_length) ; 01333 01334 /* === Kill the pivotal column ====================================== */ 01335 01336 KILL_PRINCIPAL_COL (pivot_col) ; 01337 01338 /* === Clear mark =================================================== */ 01339 01340 tag_mark += (max_deg + 1) ; 01341 if (tag_mark >= max_mark) 01342 { 01343 COLAMD_DEBUG2 (("clearing tag_mark\n")) ; 01344 tag_mark = Eigen::internal::clear_mark (n_row, Row) ; 01345 } 01346 01347 /* === Finalize the new pivot row, and column scores ================ */ 01348 01349 COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ; 01350 01351 /* for each column in pivot row */ 01352 rp = &A [pivot_row_start] ; 01353 /* compact the pivot row */ 01354 new_rp = rp ; 01355 rp_end = rp + pivot_row_length ; 01356 while (rp < rp_end) 01357 { 01358 col = *rp++ ; 01359 /* skip dead columns */ 01360 if (COL_IS_DEAD (col)) 01361 { 01362 continue ; 01363 } 01364 *new_rp++ = col ; 01365 /* add new pivot row to column */ 01366 A [Col [col].start + (Col [col].length++)] = pivot_row ; 01367 01368 /* retrieve score so far and add on pivot row's degree. */ 01369 /* (we wait until here for this in case the pivot */ 01370 /* row's degree was reduced due to mass elimination). */ 01371 cur_score = Col [col].shared2.score + pivot_row_degree ; 01372 01373 /* calculate the max possible score as the number of */ 01374 /* external columns minus the 'k' value minus the */ 01375 /* columns thickness */ 01376 max_score = n_col - k - Col [col].shared1.thickness ; 01377 01378 /* make the score the external degree of the union-of-rows */ 01379 cur_score -= Col [col].shared1.thickness ; 01380 01381 /* make sure score is less or equal than the max score */ 01382 cur_score = numext::mini(cur_score, max_score) ; 01383 COLAMD_ASSERT (cur_score >= 0) ; 01384 01385 /* store updated score */ 01386 Col [col].shared2.score = cur_score ; 01387 01388 /* === Place column back in degree list ========================= */ 01389 01390 COLAMD_ASSERT (min_score >= 0) ; 01391 COLAMD_ASSERT (min_score <= n_col) ; 01392 COLAMD_ASSERT (cur_score >= 0) ; 01393 COLAMD_ASSERT (cur_score <= n_col) ; 01394 COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ; 01395 next_col = head [cur_score] ; 01396 Col [col].shared4.degree_next = next_col ; 01397 Col [col].shared3.prev = COLAMD_EMPTY ; 01398 if (next_col != COLAMD_EMPTY) 01399 { 01400 Col [next_col].shared3.prev = col ; 01401 } 01402 head [cur_score] = col ; 01403 01404 /* see if this score is less than current min */ 01405 min_score = numext::mini(min_score, cur_score) ; 01406 01407 } 01408 01409 /* === Resurrect the new pivot row ================================== */ 01410 01411 if (pivot_row_degree > 0) 01412 { 01413 /* update pivot row length to reflect any cols that were killed */ 01414 /* during super-col detection and mass elimination */ 01415 Row [pivot_row].start = pivot_row_start ; 01416 Row [pivot_row].length = (IndexType) (new_rp - &A[pivot_row_start]) ; 01417 Row [pivot_row].shared1.degree = pivot_row_degree ; 01418 Row [pivot_row].shared2.mark = 0 ; 01419 /* pivot row is no longer dead */ 01420 } 01421 } 01422 01423 /* === All principal columns have now been ordered ====================== */ 01424 01425 return (ngarbage) ; 01426 } 01427 01428 01429 /* ========================================================================== */ 01430 /* === order_children ======================================================= */ 01431 /* ========================================================================== */ 01432 01433 /* 01434 The find_ordering routine has ordered all of the principal columns (the 01435 representatives of the supercolumns). The non-principal columns have not 01436 yet been ordered. This routine orders those columns by walking up the 01437 parent tree (a column is a child of the column which absorbed it). The 01438 final permutation vector is then placed in p [0 ... n_col-1], with p [0] 01439 being the first column, and p [n_col-1] being the last. It doesn't look 01440 like it at first glance, but be assured that this routine takes time linear 01441 in the number of columns. Although not immediately obvious, the time 01442 taken by this routine is O (n_col), that is, linear in the number of 01443 columns. Not user-callable. 01444 */ 01445 template <typename IndexType> 01446 static inline void order_children 01447 ( 01448 /* === Parameters ======================================================= */ 01449 01450 IndexType n_col, /* number of columns of A */ 01451 colamd_col<IndexType> Col [], /* of size n_col+1 */ 01452 IndexType p [] /* p [0 ... n_col-1] is the column permutation*/ 01453 ) 01454 { 01455 /* === Local variables ================================================== */ 01456 01457 IndexType i ; /* loop counter for all columns */ 01458 IndexType c ; /* column index */ 01459 IndexType parent ; /* index of column's parent */ 01460 IndexType order ; /* column's order */ 01461 01462 /* === Order each non-principal column ================================== */ 01463 01464 for (i = 0 ; i < n_col ; i++) 01465 { 01466 /* find an un-ordered non-principal column */ 01467 COLAMD_ASSERT (COL_IS_DEAD (i)) ; 01468 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY) 01469 { 01470 parent = i ; 01471 /* once found, find its principal parent */ 01472 do 01473 { 01474 parent = Col [parent].shared1.parent ; 01475 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; 01476 01477 /* now, order all un-ordered non-principal columns along path */ 01478 /* to this parent. collapse tree at the same time */ 01479 c = i ; 01480 /* get order of parent */ 01481 order = Col [parent].shared2.order ; 01482 01483 do 01484 { 01485 COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ; 01486 01487 /* order this column */ 01488 Col [c].shared2.order = order++ ; 01489 /* collaps tree */ 01490 Col [c].shared1.parent = parent ; 01491 01492 /* get immediate parent of this column */ 01493 c = Col [c].shared1.parent ; 01494 01495 /* continue until we hit an ordered column. There are */ 01496 /* guarranteed not to be anymore unordered columns */ 01497 /* above an ordered column */ 01498 } while (Col [c].shared2.order == COLAMD_EMPTY) ; 01499 01500 /* re-order the super_col parent to largest order for this group */ 01501 Col [parent].shared2.order = order ; 01502 } 01503 } 01504 01505 /* === Generate the permutation ========================================= */ 01506 01507 for (c = 0 ; c < n_col ; c++) 01508 { 01509 p [Col [c].shared2.order] = c ; 01510 } 01511 } 01512 01513 01514 /* ========================================================================== */ 01515 /* === detect_super_cols ==================================================== */ 01516 /* ========================================================================== */ 01517 01518 /* 01519 Detects supercolumns by finding matches between columns in the hash buckets. 01520 Check amongst columns in the set A [row_start ... row_start + row_length-1]. 01521 The columns under consideration are currently *not* in the degree lists, 01522 and have already been placed in the hash buckets. 01523 01524 The hash bucket for columns whose hash function is equal to h is stored 01525 as follows: 01526 01527 if head [h] is >= 0, then head [h] contains a degree list, so: 01528 01529 head [h] is the first column in degree bucket h. 01530 Col [head [h]].headhash gives the first column in hash bucket h. 01531 01532 otherwise, the degree list is empty, and: 01533 01534 -(head [h] + 2) is the first column in hash bucket h. 01535 01536 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous 01537 column" pointer. Col [c].shared3.hash is used instead as the hash number 01538 for that column. The value of Col [c].shared4.hash_next is the next column 01539 in the same hash bucket. 01540 01541 Assuming no, or "few" hash collisions, the time taken by this routine is 01542 linear in the sum of the sizes (lengths) of each column whose score has 01543 just been computed in the approximate degree computation. 01544 Not user-callable. 01545 */ 01546 template <typename IndexType> 01547 static void detect_super_cols 01548 ( 01549 /* === Parameters ======================================================= */ 01550 01551 colamd_col<IndexType> Col [], /* of size n_col+1 */ 01552 IndexType A [], /* row indices of A */ 01553 IndexType head [], /* head of degree lists and hash buckets */ 01554 IndexType row_start, /* pointer to set of columns to check */ 01555 IndexType row_length /* number of columns to check */ 01556 ) 01557 { 01558 /* === Local variables ================================================== */ 01559 01560 IndexType hash ; /* hash value for a column */ 01561 IndexType *rp ; /* pointer to a row */ 01562 IndexType c ; /* a column index */ 01563 IndexType super_c ; /* column index of the column to absorb into */ 01564 IndexType *cp1 ; /* column pointer for column super_c */ 01565 IndexType *cp2 ; /* column pointer for column c */ 01566 IndexType length ; /* length of column super_c */ 01567 IndexType prev_c ; /* column preceding c in hash bucket */ 01568 IndexType i ; /* loop counter */ 01569 IndexType *rp_end ; /* pointer to the end of the row */ 01570 IndexType col ; /* a column index in the row to check */ 01571 IndexType head_column ; /* first column in hash bucket or degree list */ 01572 IndexType first_col ; /* first column in hash bucket */ 01573 01574 /* === Consider each column in the row ================================== */ 01575 01576 rp = &A [row_start] ; 01577 rp_end = rp + row_length ; 01578 while (rp < rp_end) 01579 { 01580 col = *rp++ ; 01581 if (COL_IS_DEAD (col)) 01582 { 01583 continue ; 01584 } 01585 01586 /* get hash number for this column */ 01587 hash = Col [col].shared3.hash ; 01588 COLAMD_ASSERT (hash <= n_col) ; 01589 01590 /* === Get the first column in this hash bucket ===================== */ 01591 01592 head_column = head [hash] ; 01593 if (head_column > COLAMD_EMPTY) 01594 { 01595 first_col = Col [head_column].shared3.headhash ; 01596 } 01597 else 01598 { 01599 first_col = - (head_column + 2) ; 01600 } 01601 01602 /* === Consider each column in the hash bucket ====================== */ 01603 01604 for (super_c = first_col ; super_c != COLAMD_EMPTY ; 01605 super_c = Col [super_c].shared4.hash_next) 01606 { 01607 COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ; 01608 COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ; 01609 length = Col [super_c].length ; 01610 01611 /* prev_c is the column preceding column c in the hash bucket */ 01612 prev_c = super_c ; 01613 01614 /* === Compare super_c with all columns after it ================ */ 01615 01616 for (c = Col [super_c].shared4.hash_next ; 01617 c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next) 01618 { 01619 COLAMD_ASSERT (c != super_c) ; 01620 COLAMD_ASSERT (COL_IS_ALIVE (c)) ; 01621 COLAMD_ASSERT (Col [c].shared3.hash == hash) ; 01622 01623 /* not identical if lengths or scores are different */ 01624 if (Col [c].length != length || 01625 Col [c].shared2.score != Col [super_c].shared2.score) 01626 { 01627 prev_c = c ; 01628 continue ; 01629 } 01630 01631 /* compare the two columns */ 01632 cp1 = &A [Col [super_c].start] ; 01633 cp2 = &A [Col [c].start] ; 01634 01635 for (i = 0 ; i < length ; i++) 01636 { 01637 /* the columns are "clean" (no dead rows) */ 01638 COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ; 01639 COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ; 01640 /* row indices will same order for both supercols, */ 01641 /* no gather scatter nessasary */ 01642 if (*cp1++ != *cp2++) 01643 { 01644 break ; 01645 } 01646 } 01647 01648 /* the two columns are different if the for-loop "broke" */ 01649 if (i != length) 01650 { 01651 prev_c = c ; 01652 continue ; 01653 } 01654 01655 /* === Got it! two columns are identical =================== */ 01656 01657 COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; 01658 01659 Col [super_c].shared1.thickness += Col [c].shared1.thickness ; 01660 Col [c].shared1.parent = super_c ; 01661 KILL_NON_PRINCIPAL_COL (c) ; 01662 /* order c later, in order_children() */ 01663 Col [c].shared2.order = COLAMD_EMPTY ; 01664 /* remove c from hash bucket */ 01665 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; 01666 } 01667 } 01668 01669 /* === Empty this hash bucket ======================================= */ 01670 01671 if (head_column > COLAMD_EMPTY) 01672 { 01673 /* corresponding degree list "hash" is not empty */ 01674 Col [head_column].shared3.headhash = COLAMD_EMPTY ; 01675 } 01676 else 01677 { 01678 /* corresponding degree list "hash" is empty */ 01679 head [hash] = COLAMD_EMPTY ; 01680 } 01681 } 01682 } 01683 01684 01685 /* ========================================================================== */ 01686 /* === garbage_collection =================================================== */ 01687 /* ========================================================================== */ 01688 01689 /* 01690 Defragments and compacts columns and rows in the workspace A. Used when 01691 all avaliable memory has been used while performing row merging. Returns 01692 the index of the first free position in A, after garbage collection. The 01693 time taken by this routine is linear is the size of the array A, which is 01694 itself linear in the number of nonzeros in the input matrix. 01695 Not user-callable. 01696 */ 01697 template <typename IndexType> 01698 static IndexType garbage_collection /* returns the new value of pfree */ 01699 ( 01700 /* === Parameters ======================================================= */ 01701 01702 IndexType n_row, /* number of rows */ 01703 IndexType n_col, /* number of columns */ 01704 Colamd_Row<IndexType> Row [], /* row info */ 01705 colamd_col<IndexType> Col [], /* column info */ 01706 IndexType A [], /* A [0 ... Alen-1] holds the matrix */ 01707 IndexType *pfree /* &A [0] ... pfree is in use */ 01708 ) 01709 { 01710 /* === Local variables ================================================== */ 01711 01712 IndexType *psrc ; /* source pointer */ 01713 IndexType *pdest ; /* destination pointer */ 01714 IndexType j ; /* counter */ 01715 IndexType r ; /* a row index */ 01716 IndexType c ; /* a column index */ 01717 IndexType length ; /* length of a row or column */ 01718 01719 /* === Defragment the columns =========================================== */ 01720 01721 pdest = &A[0] ; 01722 for (c = 0 ; c < n_col ; c++) 01723 { 01724 if (COL_IS_ALIVE (c)) 01725 { 01726 psrc = &A [Col [c].start] ; 01727 01728 /* move and compact the column */ 01729 COLAMD_ASSERT (pdest <= psrc) ; 01730 Col [c].start = (IndexType) (pdest - &A [0]) ; 01731 length = Col [c].length ; 01732 for (j = 0 ; j < length ; j++) 01733 { 01734 r = *psrc++ ; 01735 if (ROW_IS_ALIVE (r)) 01736 { 01737 *pdest++ = r ; 01738 } 01739 } 01740 Col [c].length = (IndexType) (pdest - &A [Col [c].start]) ; 01741 } 01742 } 01743 01744 /* === Prepare to defragment the rows =================================== */ 01745 01746 for (r = 0 ; r < n_row ; r++) 01747 { 01748 if (ROW_IS_ALIVE (r)) 01749 { 01750 if (Row [r].length == 0) 01751 { 01752 /* this row is of zero length. cannot compact it, so kill it */ 01753 COLAMD_DEBUG3 (("Defrag row kill\n")) ; 01754 KILL_ROW (r) ; 01755 } 01756 else 01757 { 01758 /* save first column index in Row [r].shared2.first_column */ 01759 psrc = &A [Row [r].start] ; 01760 Row [r].shared2.first_column = *psrc ; 01761 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; 01762 /* flag the start of the row with the one's complement of row */ 01763 *psrc = ONES_COMPLEMENT (r) ; 01764 01765 } 01766 } 01767 } 01768 01769 /* === Defragment the rows ============================================== */ 01770 01771 psrc = pdest ; 01772 while (psrc < pfree) 01773 { 01774 /* find a negative number ... the start of a row */ 01775 if (*psrc++ < 0) 01776 { 01777 psrc-- ; 01778 /* get the row index */ 01779 r = ONES_COMPLEMENT (*psrc) ; 01780 COLAMD_ASSERT (r >= 0 && r < n_row) ; 01781 /* restore first column index */ 01782 *psrc = Row [r].shared2.first_column ; 01783 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; 01784 01785 /* move and compact the row */ 01786 COLAMD_ASSERT (pdest <= psrc) ; 01787 Row [r].start = (IndexType) (pdest - &A [0]) ; 01788 length = Row [r].length ; 01789 for (j = 0 ; j < length ; j++) 01790 { 01791 c = *psrc++ ; 01792 if (COL_IS_ALIVE (c)) 01793 { 01794 *pdest++ = c ; 01795 } 01796 } 01797 Row [r].length = (IndexType) (pdest - &A [Row [r].start]) ; 01798 01799 } 01800 } 01801 /* ensure we found all the rows */ 01802 COLAMD_ASSERT (debug_rows == 0) ; 01803 01804 /* === Return the new value of pfree ==================================== */ 01805 01806 return ((IndexType) (pdest - &A [0])) ; 01807 } 01808 01809 01810 /* ========================================================================== */ 01811 /* === clear_mark =========================================================== */ 01812 /* ========================================================================== */ 01813 01814 /* 01815 Clears the Row [].shared2.mark array, and returns the new tag_mark. 01816 Return value is the new tag_mark. Not user-callable. 01817 */ 01818 template <typename IndexType> 01819 static inline IndexType clear_mark /* return the new value for tag_mark */ 01820 ( 01821 /* === Parameters ======================================================= */ 01822 01823 IndexType n_row, /* number of rows in A */ 01824 Colamd_Row<IndexType> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ 01825 ) 01826 { 01827 /* === Local variables ================================================== */ 01828 01829 IndexType r ; 01830 01831 for (r = 0 ; r < n_row ; r++) 01832 { 01833 if (ROW_IS_ALIVE (r)) 01834 { 01835 Row [r].shared2.mark = 0 ; 01836 } 01837 } 01838 return (1) ; 01839 } 01840 01841 01842 } // namespace internal 01843 #endif