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