| // // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| // This file is modified from the colamd/symamd library. The copyright is below |
| |
| // The authors of the code itself are Stefan I. Larimore and Timothy A. |
| // Davis (davis@cise.ufl.edu), University of Florida. The algorithm was |
| // developed in collaboration with John Gilbert, Xerox PARC, and Esmond |
| // Ng, Oak Ridge National Laboratory. |
| // |
| // Date: |
| // |
| // September 8, 2003. Version 2.3. |
| // |
| // Acknowledgements: |
| // |
| // This work was supported by the National Science Foundation, under |
| // grants DMS-9504974 and DMS-9803599. |
| // |
| // Notice: |
| // |
| // Copyright (c) 1998-2003 by the University of Florida. |
| // All Rights Reserved. |
| // |
| // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
| // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
| // |
| // Permission is hereby granted to use, copy, modify, and/or distribute |
| // this program, provided that the Copyright, this License, and the |
| // Availability of the original version is retained on all copies and made |
| // accessible to the end-user of any code or package that includes COLAMD |
| // or any modified version of COLAMD. |
| // |
| // Availability: |
| // |
| // The colamd/symamd library is available at |
| // |
| // http://www.suitesparse.com |
| |
| #ifndef EIGEN_COLAMD_H |
| #define EIGEN_COLAMD_H |
| |
| namespace internal { |
| |
| namespace Colamd { |
| |
| /* Ensure that debugging is turned off: */ |
| #ifndef COLAMD_NDEBUG |
| #define COLAMD_NDEBUG |
| #endif /* NDEBUG */ |
| |
| /* ========================================================================== */ |
| /* === Knob and statistics definitions ====================================== */ |
| /* ========================================================================== */ |
| |
| /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ |
| const int NKnobs = 20; |
| |
| /* number of output statistics. Only stats [0..6] are currently used. */ |
| const int NStats = 20; |
| |
| /* Indices into knobs and stats array. */ |
| enum KnobsStatsIndex { |
| /* knobs [0] and stats [0]: dense row knob and output statistic. */ |
| DenseRow = 0, |
| |
| /* knobs [1] and stats [1]: dense column knob and output statistic. */ |
| DenseCol = 1, |
| |
| /* stats [2]: memory defragmentation count output statistic */ |
| DefragCount = 2, |
| |
| /* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ |
| Status = 3, |
| |
| /* stats [4..6]: error info, or info on jumbled columns */ |
| Info1 = 4, |
| Info2 = 5, |
| Info3 = 6 |
| }; |
| |
| /* error codes returned in stats [3]: */ |
| enum Status { |
| Ok = 0, |
| OkButJumbled = 1, |
| ErrorANotPresent = -1, |
| ErrorPNotPresent = -2, |
| ErrorNrowNegative = -3, |
| ErrorNcolNegative = -4, |
| ErrorNnzNegative = -5, |
| ErrorP0Nonzero = -6, |
| ErrorATooSmall = -7, |
| ErrorColLengthNegative = -8, |
| ErrorRowIndexOutOfBounds = -9, |
| ErrorOutOfMemory = -10, |
| ErrorInternalError = -999 |
| }; |
| /* ========================================================================== */ |
| /* === Definitions ========================================================== */ |
| /* ========================================================================== */ |
| |
| template <typename IndexType> |
| IndexType ones_complement(const IndexType r) { |
| return (-(r)-1); |
| } |
| |
| /* -------------------------------------------------------------------------- */ |
| const int Empty = -1; |
| |
| /* Row and column status */ |
| enum RowColumnStatus { Alive = 0, Dead = -1 }; |
| |
| /* Column status */ |
| enum ColumnStatus { DeadPrincipal = -1, DeadNonPrincipal = -2 }; |
| |
| /* ========================================================================== */ |
| /* === Colamd reporting mechanism =========================================== */ |
| /* ========================================================================== */ |
| |
| // == Row and Column structures == |
| template <typename IndexType> |
| struct ColStructure { |
| IndexType start; /* index for A of first row in this column, or Dead */ |
| /* if column is dead */ |
| IndexType length; /* number of rows in this column */ |
| union { |
| IndexType thickness; /* number of original columns represented by this */ |
| /* col, if the column is alive */ |
| IndexType parent; /* parent in parent tree super-column structure, if */ |
| /* the column is dead */ |
| } shared1; |
| union { |
| IndexType score; /* the score used to maintain heap, if col is alive */ |
| IndexType order; /* pivot ordering of this column, if col is dead */ |
| } shared2; |
| union { |
| IndexType headhash; /* head of a hash bucket, if col is at the head of */ |
| /* a degree list */ |
| IndexType hash; /* hash value, if col is not in a degree list */ |
| IndexType prev; /* previous column in degree list, if col is in a */ |
| /* degree list (but not at the head of a degree list) */ |
| } shared3; |
| union { |
| IndexType degree_next; /* next column, if col is in a degree list */ |
| IndexType hash_next; /* next column, if col is in a hash list */ |
| } shared4; |
| |
| inline bool is_dead() const { return start < Alive; } |
| |
| inline bool is_alive() const { return start >= Alive; } |
| |
| inline bool is_dead_principal() const { return start == DeadPrincipal; } |
| |
| inline void kill_principal() { start = DeadPrincipal; } |
| |
| inline void kill_non_principal() { start = DeadNonPrincipal; } |
| }; |
| |
| template <typename IndexType> |
| struct RowStructure { |
| IndexType start; /* index for A of first col in this row */ |
| IndexType length; /* number of principal columns in this row */ |
| union { |
| IndexType degree; /* number of principal & non-principal columns in row */ |
| IndexType p; /* used as a row pointer in init_rows_cols () */ |
| } shared1; |
| union { |
| IndexType mark; /* for computing set differences and marking dead rows*/ |
| IndexType first_column; /* first column in row (used in garbage collection) */ |
| } shared2; |
| |
| inline bool is_dead() const { return shared2.mark < Alive; } |
| |
| inline bool is_alive() const { return shared2.mark >= Alive; } |
| |
| inline void kill() { shared2.mark = Dead; } |
| }; |
| |
| /* ========================================================================== */ |
| /* === Colamd recommended memory size ======================================= */ |
| /* ========================================================================== */ |
| |
| /* |
| The recommended length Alen of the array A passed to colamd is given by |
| the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any |
| argument is negative. 2*nnz space is required for the row and column |
| indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is |
| required for the Col and Row arrays, respectively, which are internal to |
| colamd. An additional n_col space is the minimal amount of "elbow room", |
| and nnz/5 more space is recommended for run time efficiency. |
| |
| This macro is not needed when using symamd. |
| |
| Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid |
| gcc -pedantic warning messages. |
| */ |
| template <typename IndexType> |
| inline IndexType colamd_c(IndexType n_col) { |
| return IndexType(((n_col) + 1) * sizeof(ColStructure<IndexType>) / sizeof(IndexType)); |
| } |
| |
| template <typename IndexType> |
| inline IndexType colamd_r(IndexType n_row) { |
| return IndexType(((n_row) + 1) * sizeof(RowStructure<IndexType>) / sizeof(IndexType)); |
| } |
| |
| // Prototypes of non-user callable routines |
| template <typename IndexType> |
| static IndexType init_rows_cols(IndexType n_row, IndexType n_col, RowStructure<IndexType> Row[], |
| ColStructure<IndexType> col[], IndexType A[], IndexType p[], IndexType stats[NStats]); |
| |
| template <typename IndexType> |
| static void init_scoring(IndexType n_row, IndexType n_col, RowStructure<IndexType> Row[], ColStructure<IndexType> Col[], |
| IndexType A[], IndexType head[], double knobs[NKnobs], IndexType *p_n_row2, |
| IndexType *p_n_col2, IndexType *p_max_deg); |
| |
| template <typename IndexType> |
| static IndexType find_ordering(IndexType n_row, IndexType n_col, IndexType Alen, RowStructure<IndexType> Row[], |
| ColStructure<IndexType> Col[], IndexType A[], IndexType head[], IndexType n_col2, |
| IndexType max_deg, IndexType pfree); |
| |
| template <typename IndexType> |
| static void order_children(IndexType n_col, ColStructure<IndexType> Col[], IndexType p[]); |
| |
| template <typename IndexType> |
| static void detect_super_cols(ColStructure<IndexType> Col[], IndexType A[], IndexType head[], IndexType row_start, |
| IndexType row_length); |
| |
| template <typename IndexType> |
| static IndexType garbage_collection(IndexType n_row, IndexType n_col, RowStructure<IndexType> Row[], |
| ColStructure<IndexType> Col[], IndexType A[], IndexType *pfree); |
| |
| template <typename IndexType> |
| static inline IndexType clear_mark(IndexType n_row, RowStructure<IndexType> Row[]); |
| |
| /* === No debugging ========================================================= */ |
| |
| #define COLAMD_DEBUG0(params) ; |
| #define COLAMD_DEBUG1(params) ; |
| #define COLAMD_DEBUG2(params) ; |
| #define COLAMD_DEBUG3(params) ; |
| #define COLAMD_DEBUG4(params) ; |
| |
| #define COLAMD_ASSERT(expression) ((void)0) |
| |
| /** |
| * \brief Returns the recommended value of Alen |
| * |
| * Returns recommended value of Alen for use by colamd. |
| * Returns -1 if any input argument is negative. |
| * The use of this routine or macro is optional. |
| * Note that the macro uses its arguments more than once, |
| * so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED. |
| * |
| * \param nnz nonzeros in A |
| * \param n_row number of rows in A |
| * \param n_col number of columns in A |
| * \return recommended value of Alen for use by colamd |
| */ |
| template <typename IndexType> |
| inline IndexType recommended(IndexType nnz, IndexType n_row, IndexType n_col) { |
| if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) |
| return (-1); |
| else |
| return (2 * (nnz) + colamd_c(n_col) + colamd_r(n_row) + (n_col) + ((nnz) / 5)); |
| } |
| |
| /** |
| * \brief set default parameters The use of this routine is optional. |
| * |
| * Colamd: rows with more than (knobs [DenseRow] * n_col) |
| * entries are removed prior to ordering. Columns with more than |
| * (knobs [DenseCol] * n_row) entries are removed prior to |
| * ordering, and placed last in the output column ordering. |
| * |
| * DenseRow and DenseCol are defined as 0 and 1, |
| * respectively, in colamd.h. Default values of these two knobs |
| * are both 0.5. Currently, only knobs [0] and knobs [1] are |
| * used, but future versions may use more knobs. If so, they will |
| * be properly set to their defaults by the future version of |
| * colamd_set_defaults, so that the code that calls colamd will |
| * not need to change, assuming that you either use |
| * colamd_set_defaults, or pass a (double *) NULL pointer as the |
| * knobs array to colamd or symamd. |
| * |
| * \param knobs parameter settings for colamd |
| */ |
| |
| static inline void set_defaults(double knobs[NKnobs]) { |
| /* === Local variables ================================================== */ |
| |
| int i; |
| |
| if (!knobs) { |
| return; /* no knobs to initialize */ |
| } |
| for (i = 0; i < NKnobs; i++) { |
| knobs[i] = 0; |
| } |
| knobs[Colamd::DenseRow] = 0.5; /* ignore rows over 50% dense */ |
| knobs[Colamd::DenseCol] = 0.5; /* ignore columns over 50% dense */ |
| } |
| |
| /** |
| * \brief Computes a column ordering using the column approximate minimum degree ordering |
| * |
| * Computes a column ordering (Q) of A such that P(AQ)=LU or |
| * (AQ)'AQ=LL' have less fill-in and require fewer floating point |
| * operations than factorizing the unpermuted matrix A or A'A, |
| * respectively. |
| * |
| * |
| * \param n_row number of rows in A |
| * \param n_col number of columns in A |
| * \param Alen, size of the array A |
| * \param A row indices of the matrix, of size ALen |
| * \param p column pointers of A, of size n_col+1 |
| * \param knobs parameter settings for colamd |
| * \param stats colamd output statistics and error codes |
| */ |
| template <typename IndexType> |
| static bool compute_ordering(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, |
| double knobs[NKnobs], IndexType stats[NStats]) { |
| /* === Local variables ================================================== */ |
| |
| IndexType i; /* loop index */ |
| IndexType nnz; /* nonzeros in A */ |
| IndexType Row_size; /* size of Row [], in integers */ |
| IndexType Col_size; /* size of Col [], in integers */ |
| IndexType need; /* minimum required length of A */ |
| Colamd::RowStructure<IndexType> *Row; /* pointer into A of Row [0..n_row] array */ |
| Colamd::ColStructure<IndexType> *Col; /* pointer into A of Col [0..n_col] array */ |
| IndexType n_col2; /* number of non-dense, non-empty columns */ |
| IndexType n_row2; /* number of non-dense, non-empty rows */ |
| IndexType ngarbage; /* number of garbage collections performed */ |
| IndexType max_deg; /* maximum row degree */ |
| double default_knobs[NKnobs]; /* default knobs array */ |
| |
| /* === Check the input arguments ======================================== */ |
| |
| if (!stats) { |
| COLAMD_DEBUG0(("colamd: stats not present\n")); |
| return (false); |
| } |
| for (i = 0; i < NStats; i++) { |
| stats[i] = 0; |
| } |
| stats[Colamd::Status] = Colamd::Ok; |
| stats[Colamd::Info1] = -1; |
| stats[Colamd::Info2] = -1; |
| |
| if (!A) /* A is not present */ |
| { |
| stats[Colamd::Status] = Colamd::ErrorANotPresent; |
| COLAMD_DEBUG0(("colamd: A not present\n")); |
| return (false); |
| } |
| |
| if (!p) /* p is not present */ |
| { |
| stats[Colamd::Status] = Colamd::ErrorPNotPresent; |
| COLAMD_DEBUG0(("colamd: p not present\n")); |
| return (false); |
| } |
| |
| if (n_row < 0) /* n_row must be >= 0 */ |
| { |
| stats[Colamd::Status] = Colamd::ErrorNrowNegative; |
| stats[Colamd::Info1] = n_row; |
| COLAMD_DEBUG0(("colamd: nrow negative %d\n", n_row)); |
| return (false); |
| } |
| |
| if (n_col < 0) /* n_col must be >= 0 */ |
| { |
| stats[Colamd::Status] = Colamd::ErrorNcolNegative; |
| stats[Colamd::Info1] = n_col; |
| COLAMD_DEBUG0(("colamd: ncol negative %d\n", n_col)); |
| return (false); |
| } |
| |
| nnz = p[n_col]; |
| if (nnz < 0) /* nnz must be >= 0 */ |
| { |
| stats[Colamd::Status] = Colamd::ErrorNnzNegative; |
| stats[Colamd::Info1] = nnz; |
| COLAMD_DEBUG0(("colamd: number of entries negative %d\n", nnz)); |
| return (false); |
| } |
| |
| if (p[0] != 0) { |
| stats[Colamd::Status] = Colamd::ErrorP0Nonzero; |
| stats[Colamd::Info1] = p[0]; |
| COLAMD_DEBUG0(("colamd: p[0] not zero %d\n", p[0])); |
| return (false); |
| } |
| |
| /* === If no knobs, set default knobs =================================== */ |
| |
| if (!knobs) { |
| set_defaults(default_knobs); |
| knobs = default_knobs; |
| } |
| |
| /* === Allocate the Row and Col arrays from array A ===================== */ |
| |
| Col_size = colamd_c(n_col); |
| Row_size = colamd_r(n_row); |
| need = 2 * nnz + n_col + Col_size + Row_size; |
| |
| if (need > Alen) { |
| /* not enough space in array A to perform the ordering */ |
| stats[Colamd::Status] = Colamd::ErrorATooSmall; |
| stats[Colamd::Info1] = need; |
| stats[Colamd::Info2] = Alen; |
| COLAMD_DEBUG0(("colamd: Need Alen >= %d, given only Alen = %d\n", need, Alen)); |
| return (false); |
| } |
| |
| Alen -= Col_size + Row_size; |
| Col = (ColStructure<IndexType> *)&A[Alen]; |
| Row = (RowStructure<IndexType> *)&A[Alen + Col_size]; |
| |
| /* === Construct the row and column data structures ===================== */ |
| |
| if (!Colamd::init_rows_cols(n_row, n_col, Row, Col, A, p, stats)) { |
| /* input matrix is invalid */ |
| COLAMD_DEBUG0(("colamd: Matrix invalid\n")); |
| return (false); |
| } |
| |
| /* === Initialize scores, kill dense rows/columns ======================= */ |
| |
| Colamd::init_scoring(n_row, n_col, Row, Col, A, p, knobs, &n_row2, &n_col2, &max_deg); |
| |
| /* === Order the supercolumns =========================================== */ |
| |
| ngarbage = Colamd::find_ordering(n_row, n_col, Alen, Row, Col, A, p, n_col2, max_deg, 2 * nnz); |
| |
| /* === Order the non-principal columns ================================== */ |
| |
| Colamd::order_children(n_col, Col, p); |
| |
| /* === Return statistics in stats ======================================= */ |
| |
| stats[Colamd::DenseRow] = n_row - n_row2; |
| stats[Colamd::DenseCol] = n_col - n_col2; |
| stats[Colamd::DefragCount] = ngarbage; |
| COLAMD_DEBUG0(("colamd: done.\n")); |
| return (true); |
| } |
| |
| /* ========================================================================== */ |
| /* === NON-USER-CALLABLE ROUTINES: ========================================== */ |
| /* ========================================================================== */ |
| |
| /* There are no user-callable routines beyond this point in the file */ |
| |
| /* ========================================================================== */ |
| /* === init_rows_cols ======================================================= */ |
| /* ========================================================================== */ |
| |
| /* |
| Takes the column form of the matrix in A and creates the row form of the |
| matrix. Also, row and column attributes are stored in the Col and Row |
| structs. If the columns are un-sorted or contain duplicate row indices, |
| this routine will also sort and remove duplicate row indices from the |
| column form of the matrix. Returns false if the matrix is invalid, |
| true otherwise. Not user-callable. |
| */ |
| template <typename IndexType> |
| static IndexType init_rows_cols /* returns true if OK, or false otherwise */ |
| ( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_row, /* number of rows of A */ |
| IndexType n_col, /* number of columns of A */ |
| RowStructure<IndexType> Row[], /* of size n_row+1 */ |
| ColStructure<IndexType> Col[], /* of size n_col+1 */ |
| IndexType A[], /* row indices of A, of size Alen */ |
| IndexType p[], /* pointers to columns in A, of size n_col+1 */ |
| IndexType stats[NStats] /* colamd statistics */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType col; /* a column index */ |
| IndexType row; /* a row index */ |
| IndexType *cp; /* a column pointer */ |
| IndexType *cp_end; /* a pointer to the end of a column */ |
| IndexType *rp; /* a row pointer */ |
| IndexType *rp_end; /* a pointer to the end of a row */ |
| IndexType last_row; /* previous row */ |
| |
| /* === Initialize columns, and check column pointers ==================== */ |
| |
| for (col = 0; col < n_col; col++) { |
| Col[col].start = p[col]; |
| Col[col].length = p[col + 1] - p[col]; |
| |
| if ((Col[col].length) < 0) // extra parentheses to work-around gcc bug 10200 |
| { |
| /* column pointers must be non-decreasing */ |
| stats[Colamd::Status] = Colamd::ErrorColLengthNegative; |
| stats[Colamd::Info1] = col; |
| stats[Colamd::Info2] = Col[col].length; |
| COLAMD_DEBUG0(("colamd: col %d length %d < 0\n", col, Col[col].length)); |
| return (false); |
| } |
| |
| Col[col].shared1.thickness = 1; |
| Col[col].shared2.score = 0; |
| Col[col].shared3.prev = Empty; |
| Col[col].shared4.degree_next = Empty; |
| } |
| |
| /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ |
| |
| /* === Scan columns, compute row degrees, and check row indices ========= */ |
| |
| stats[Info3] = 0; /* number of duplicate or unsorted row indices*/ |
| |
| for (row = 0; row < n_row; row++) { |
| Row[row].length = 0; |
| Row[row].shared2.mark = -1; |
| } |
| |
| for (col = 0; col < n_col; col++) { |
| last_row = -1; |
| |
| cp = &A[p[col]]; |
| cp_end = &A[p[col + 1]]; |
| |
| while (cp < cp_end) { |
| row = *cp++; |
| |
| /* make sure row indices within range */ |
| if (row < 0 || row >= n_row) { |
| stats[Colamd::Status] = Colamd::ErrorRowIndexOutOfBounds; |
| stats[Colamd::Info1] = col; |
| stats[Colamd::Info2] = row; |
| stats[Colamd::Info3] = n_row; |
| COLAMD_DEBUG0(("colamd: row %d col %d out of bounds\n", row, col)); |
| return (false); |
| } |
| |
| if (row <= last_row || Row[row].shared2.mark == col) { |
| /* row index are unsorted or repeated (or both), thus col */ |
| /* is jumbled. This is a notice, not an error condition. */ |
| stats[Colamd::Status] = Colamd::OkButJumbled; |
| stats[Colamd::Info1] = col; |
| stats[Colamd::Info2] = row; |
| (stats[Colamd::Info3])++; |
| COLAMD_DEBUG1(("colamd: row %d col %d unsorted/duplicate\n", row, col)); |
| } |
| |
| if (Row[row].shared2.mark != col) { |
| Row[row].length++; |
| } else { |
| /* this is a repeated entry in the column, */ |
| /* it will be removed */ |
| Col[col].length--; |
| } |
| |
| /* mark the row as having been seen in this column */ |
| Row[row].shared2.mark = col; |
| |
| last_row = row; |
| } |
| } |
| |
| /* === Compute row pointers ============================================= */ |
| |
| /* row form of the matrix starts directly after the column */ |
| /* form of matrix in A */ |
| Row[0].start = p[n_col]; |
| Row[0].shared1.p = Row[0].start; |
| Row[0].shared2.mark = -1; |
| for (row = 1; row < n_row; row++) { |
| Row[row].start = Row[row - 1].start + Row[row - 1].length; |
| Row[row].shared1.p = Row[row].start; |
| Row[row].shared2.mark = -1; |
| } |
| |
| /* === Create row form ================================================== */ |
| |
| if (stats[Status] == OkButJumbled) { |
| /* if cols jumbled, watch for repeated row indices */ |
| for (col = 0; col < n_col; col++) { |
| cp = &A[p[col]]; |
| cp_end = &A[p[col + 1]]; |
| while (cp < cp_end) { |
| row = *cp++; |
| if (Row[row].shared2.mark != col) { |
| A[(Row[row].shared1.p)++] = col; |
| Row[row].shared2.mark = col; |
| } |
| } |
| } |
| } else { |
| /* if cols not jumbled, we don't need the mark (this is faster) */ |
| for (col = 0; col < n_col; col++) { |
| cp = &A[p[col]]; |
| cp_end = &A[p[col + 1]]; |
| while (cp < cp_end) { |
| A[(Row[*cp++].shared1.p)++] = col; |
| } |
| } |
| } |
| |
| /* === Clear the row marks and set row degrees ========================== */ |
| |
| for (row = 0; row < n_row; row++) { |
| Row[row].shared2.mark = 0; |
| Row[row].shared1.degree = Row[row].length; |
| } |
| |
| /* === See if we need to re-create columns ============================== */ |
| |
| if (stats[Status] == OkButJumbled) { |
| COLAMD_DEBUG0(("colamd: reconstructing column form, matrix jumbled\n")); |
| |
| /* === Compute col pointers ========================================= */ |
| |
| /* col form of the matrix starts at A [0]. */ |
| /* Note, we may have a gap between the col form and the row */ |
| /* form if there were duplicate entries, if so, it will be */ |
| /* removed upon the first garbage collection */ |
| Col[0].start = 0; |
| p[0] = Col[0].start; |
| for (col = 1; col < n_col; col++) { |
| /* note that the lengths here are for pruned columns, i.e. */ |
| /* no duplicate row indices will exist for these columns */ |
| Col[col].start = Col[col - 1].start + Col[col - 1].length; |
| p[col] = Col[col].start; |
| } |
| |
| /* === Re-create col form =========================================== */ |
| |
| for (row = 0; row < n_row; row++) { |
| rp = &A[Row[row].start]; |
| rp_end = rp + Row[row].length; |
| while (rp < rp_end) { |
| A[(p[*rp++])++] = row; |
| } |
| } |
| } |
| |
| /* === Done. Matrix is not (or no longer) jumbled ====================== */ |
| |
| return (true); |
| } |
| |
| /* ========================================================================== */ |
| /* === init_scoring ========================================================= */ |
| /* ========================================================================== */ |
| |
| /* |
| Kills dense or empty columns and rows, calculates an initial score for |
| each column, and places all columns in the degree lists. Not user-callable. |
| */ |
| template <typename IndexType> |
| static void init_scoring( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_row, /* number of rows of A */ |
| IndexType n_col, /* number of columns of A */ |
| RowStructure<IndexType> Row[], /* of size n_row+1 */ |
| ColStructure<IndexType> Col[], /* of size n_col+1 */ |
| IndexType A[], /* column form and row form of A */ |
| IndexType head[], /* of size n_col+1 */ |
| double knobs[NKnobs], /* parameters */ |
| IndexType *p_n_row2, /* number of non-dense, non-empty rows */ |
| IndexType *p_n_col2, /* number of non-dense, non-empty columns */ |
| IndexType *p_max_deg /* maximum row degree */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType c; /* a column index */ |
| IndexType r, row; /* a row index */ |
| IndexType *cp; /* a column pointer */ |
| IndexType deg; /* degree of a row or column */ |
| IndexType *cp_end; /* a pointer to the end of a column */ |
| IndexType *new_cp; /* new column pointer */ |
| IndexType col_length; /* length of pruned column */ |
| IndexType score; /* current column score */ |
| IndexType n_col2; /* number of non-dense, non-empty columns */ |
| IndexType n_row2; /* number of non-dense, non-empty rows */ |
| IndexType dense_row_count; /* remove rows with more entries than this */ |
| IndexType dense_col_count; /* remove cols with more entries than this */ |
| IndexType min_score; /* smallest column score */ |
| IndexType max_deg; /* maximum row degree */ |
| IndexType next_col; /* Used to add to degree list.*/ |
| |
| /* === Extract knobs ==================================================== */ |
| |
| dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs[Colamd::DenseRow] * n_col), n_col)); |
| dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs[Colamd::DenseCol] * n_row), n_row)); |
| COLAMD_DEBUG1(("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)); |
| max_deg = 0; |
| n_col2 = n_col; |
| n_row2 = n_row; |
| |
| /* === Kill empty columns =============================================== */ |
| |
| /* Put the empty columns at the end in their natural order, so that LU */ |
| /* factorization can proceed as far as possible. */ |
| for (c = n_col - 1; c >= 0; c--) { |
| deg = Col[c].length; |
| if (deg == 0) { |
| /* this is a empty column, kill and order it last */ |
| Col[c].shared2.order = --n_col2; |
| Col[c].kill_principal(); |
| } |
| } |
| COLAMD_DEBUG1(("colamd: null columns killed: %d\n", n_col - n_col2)); |
| |
| /* === Kill dense columns =============================================== */ |
| |
| /* Put the dense columns at the end, in their natural order */ |
| for (c = n_col - 1; c >= 0; c--) { |
| /* skip any dead columns */ |
| if (Col[c].is_dead()) { |
| continue; |
| } |
| deg = Col[c].length; |
| if (deg > dense_col_count) { |
| /* this is a dense column, kill and order it last */ |
| Col[c].shared2.order = --n_col2; |
| /* decrement the row degrees */ |
| cp = &A[Col[c].start]; |
| cp_end = cp + Col[c].length; |
| while (cp < cp_end) { |
| Row[*cp++].shared1.degree--; |
| } |
| Col[c].kill_principal(); |
| } |
| } |
| COLAMD_DEBUG1(("colamd: Dense and null columns killed: %d\n", n_col - n_col2)); |
| |
| /* === Kill dense and empty rows ======================================== */ |
| |
| for (r = 0; r < n_row; r++) { |
| deg = Row[r].shared1.degree; |
| COLAMD_ASSERT(deg >= 0 && deg <= n_col); |
| if (deg > dense_row_count || deg == 0) { |
| /* kill a dense or empty row */ |
| Row[r].kill(); |
| --n_row2; |
| } else { |
| /* keep track of max degree of remaining rows */ |
| max_deg = numext::maxi(max_deg, deg); |
| } |
| } |
| COLAMD_DEBUG1(("colamd: Dense and null rows killed: %d\n", n_row - n_row2)); |
| |
| /* === Compute initial column scores ==================================== */ |
| |
| /* At this point the row degrees are accurate. They reflect the number */ |
| /* of "live" (non-dense) columns in each row. No empty rows exist. */ |
| /* Some "live" columns may contain only dead rows, however. These are */ |
| /* pruned in the code below. */ |
| |
| /* now find the initial matlab score for each column */ |
| for (c = n_col - 1; c >= 0; c--) { |
| /* skip dead column */ |
| if (Col[c].is_dead()) { |
| continue; |
| } |
| score = 0; |
| cp = &A[Col[c].start]; |
| new_cp = cp; |
| cp_end = cp + Col[c].length; |
| while (cp < cp_end) { |
| /* get a row */ |
| row = *cp++; |
| /* skip if dead */ |
| if (Row[row].is_dead()) { |
| continue; |
| } |
| /* compact the column */ |
| *new_cp++ = row; |
| /* add row's external degree */ |
| score += Row[row].shared1.degree - 1; |
| /* guard against integer overflow */ |
| score = numext::mini(score, n_col); |
| } |
| /* determine pruned column length */ |
| col_length = (IndexType)(new_cp - &A[Col[c].start]); |
| if (col_length == 0) { |
| /* a newly-made null column (all rows in this col are "dense" */ |
| /* and have already been killed) */ |
| COLAMD_DEBUG2(("Newly null killed: %d\n", c)); |
| Col[c].shared2.order = --n_col2; |
| Col[c].kill_principal(); |
| } else { |
| /* set column length and set score */ |
| COLAMD_ASSERT(score >= 0); |
| COLAMD_ASSERT(score <= n_col); |
| Col[c].length = col_length; |
| Col[c].shared2.score = score; |
| } |
| } |
| COLAMD_DEBUG1(("colamd: Dense, null, and newly-null columns killed: %d\n", n_col - n_col2)); |
| |
| /* At this point, all empty rows and columns are dead. All live columns */ |
| /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ |
| /* yet). Rows may contain dead columns, but all live rows contain at */ |
| /* least one live column. */ |
| |
| /* === Initialize degree lists ========================================== */ |
| |
| /* clear the hash buckets */ |
| for (c = 0; c <= n_col; c++) { |
| head[c] = Empty; |
| } |
| min_score = n_col; |
| /* place in reverse order, so low column indices are at the front */ |
| /* of the lists. This is to encourage natural tie-breaking */ |
| for (c = n_col - 1; c >= 0; c--) { |
| /* only add principal columns to degree lists */ |
| if (Col[c].is_alive()) { |
| COLAMD_DEBUG4(("place %d score %d minscore %d ncol %d\n", c, Col[c].shared2.score, min_score, n_col)); |
| |
| /* === Add columns score to DList =============================== */ |
| |
| score = Col[c].shared2.score; |
| |
| COLAMD_ASSERT(min_score >= 0); |
| COLAMD_ASSERT(min_score <= n_col); |
| COLAMD_ASSERT(score >= 0); |
| COLAMD_ASSERT(score <= n_col); |
| COLAMD_ASSERT(head[score] >= Empty); |
| |
| /* now add this column to dList at proper score location */ |
| next_col = head[score]; |
| Col[c].shared3.prev = Empty; |
| Col[c].shared4.degree_next = next_col; |
| |
| /* if there already was a column with the same score, set its */ |
| /* previous pointer to this new column */ |
| if (next_col != Empty) { |
| Col[next_col].shared3.prev = c; |
| } |
| head[score] = c; |
| |
| /* see if this score is less than current min */ |
| min_score = numext::mini(min_score, score); |
| } |
| } |
| |
| /* === Return number of remaining columns, and max row degree =========== */ |
| |
| *p_n_col2 = n_col2; |
| *p_n_row2 = n_row2; |
| *p_max_deg = max_deg; |
| } |
| |
| /* ========================================================================== */ |
| /* === find_ordering ======================================================== */ |
| /* ========================================================================== */ |
| |
| /* |
| Order the principal columns of the supercolumn form of the matrix |
| (no supercolumns on input). Uses a minimum approximate column minimum |
| degree ordering method. Not user-callable. |
| */ |
| template <typename IndexType> |
| static IndexType find_ordering /* return the number of garbage collections */ |
| ( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_row, /* number of rows of A */ |
| IndexType n_col, /* number of columns of A */ |
| IndexType Alen, /* size of A, 2*nnz + n_col or larger */ |
| RowStructure<IndexType> Row[], /* of size n_row+1 */ |
| ColStructure<IndexType> Col[], /* of size n_col+1 */ |
| IndexType A[], /* column form and row form of A */ |
| IndexType head[], /* of size n_col+1 */ |
| IndexType n_col2, /* Remaining columns to order */ |
| IndexType max_deg, /* Maximum row degree */ |
| IndexType pfree /* index of first free slot (2*nnz on entry) */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType k; /* current pivot ordering step */ |
| IndexType pivot_col; /* current pivot column */ |
| IndexType *cp; /* a column pointer */ |
| IndexType *rp; /* a row pointer */ |
| IndexType pivot_row; /* current pivot row */ |
| IndexType *new_cp; /* modified column pointer */ |
| IndexType *new_rp; /* modified row pointer */ |
| IndexType pivot_row_start; /* pointer to start of pivot row */ |
| IndexType pivot_row_degree; /* number of columns in pivot row */ |
| IndexType pivot_row_length; /* number of supercolumns in pivot row */ |
| IndexType pivot_col_score; /* score of pivot column */ |
| IndexType needed_memory; /* free space needed for pivot row */ |
| IndexType *cp_end; /* pointer to the end of a column */ |
| IndexType *rp_end; /* pointer to the end of a row */ |
| IndexType row; /* a row index */ |
| IndexType col; /* a column index */ |
| IndexType max_score; /* maximum possible score */ |
| IndexType cur_score; /* score of current column */ |
| unsigned int hash; /* hash value for supernode detection */ |
| IndexType head_column; /* head of hash bucket */ |
| IndexType first_col; /* first column in hash bucket */ |
| IndexType tag_mark; /* marker value for mark array */ |
| IndexType row_mark; /* Row [row].shared2.mark */ |
| IndexType set_difference; /* set difference size of row with pivot row */ |
| IndexType min_score; /* smallest column score */ |
| IndexType col_thickness; /* "thickness" (no. of columns in a supercol) */ |
| IndexType max_mark; /* maximum value of tag_mark */ |
| IndexType pivot_col_thickness; /* number of columns represented by pivot col */ |
| IndexType prev_col; /* Used by Dlist operations. */ |
| IndexType next_col; /* Used by Dlist operations. */ |
| IndexType ngarbage; /* number of garbage collections performed */ |
| |
| /* === Initialization and clear mark ==================================== */ |
| |
| max_mark = INT_MAX - n_col; /* INT_MAX defined in <limits.h> */ |
| tag_mark = Colamd::clear_mark(n_row, Row); |
| min_score = 0; |
| ngarbage = 0; |
| COLAMD_DEBUG1(("colamd: Ordering, n_col2=%d\n", n_col2)); |
| |
| /* === Order the columns ================================================ */ |
| |
| for (k = 0; k < n_col2; /* 'k' is incremented below */) { |
| /* === Select pivot column, and order it ============================ */ |
| |
| /* make sure degree list isn't empty */ |
| COLAMD_ASSERT(min_score >= 0); |
| COLAMD_ASSERT(min_score <= n_col); |
| COLAMD_ASSERT(head[min_score] >= Empty); |
| |
| /* get pivot column from head of minimum degree list */ |
| while (min_score < n_col && head[min_score] == Empty) { |
| min_score++; |
| } |
| pivot_col = head[min_score]; |
| COLAMD_ASSERT(pivot_col >= 0 && pivot_col <= n_col); |
| next_col = Col[pivot_col].shared4.degree_next; |
| head[min_score] = next_col; |
| if (next_col != Empty) { |
| Col[next_col].shared3.prev = Empty; |
| } |
| |
| COLAMD_ASSERT(Col[pivot_col].is_alive()); |
| COLAMD_DEBUG3(("Pivot col: %d\n", pivot_col)); |
| |
| /* remember score for defrag check */ |
| pivot_col_score = Col[pivot_col].shared2.score; |
| |
| /* the pivot column is the kth column in the pivot order */ |
| Col[pivot_col].shared2.order = k; |
| |
| /* increment order count by column thickness */ |
| pivot_col_thickness = Col[pivot_col].shared1.thickness; |
| k += pivot_col_thickness; |
| COLAMD_ASSERT(pivot_col_thickness > 0); |
| |
| /* === Garbage_collection, if necessary ============================= */ |
| |
| needed_memory = numext::mini(pivot_col_score, n_col - k); |
| if (pfree + needed_memory >= Alen) { |
| pfree = Colamd::garbage_collection(n_row, n_col, Row, Col, A, &A[pfree]); |
| ngarbage++; |
| /* after garbage collection we will have enough */ |
| COLAMD_ASSERT(pfree + needed_memory < Alen); |
| /* garbage collection has wiped out the Row[].shared2.mark array */ |
| tag_mark = Colamd::clear_mark(n_row, Row); |
| } |
| |
| /* === Compute pivot row pattern ==================================== */ |
| |
| /* get starting location for this new merged row */ |
| pivot_row_start = pfree; |
| |
| /* initialize new row counts to zero */ |
| pivot_row_degree = 0; |
| |
| /* tag pivot column as having been visited so it isn't included */ |
| /* in merged pivot row */ |
| Col[pivot_col].shared1.thickness = -pivot_col_thickness; |
| |
| /* pivot row is the union of all rows in the pivot column pattern */ |
| cp = &A[Col[pivot_col].start]; |
| cp_end = cp + Col[pivot_col].length; |
| while (cp < cp_end) { |
| /* get a row */ |
| row = *cp++; |
| COLAMD_DEBUG4(("Pivot col pattern %d %d\n", Row[row].is_alive(), row)); |
| /* skip if row is dead */ |
| if (Row[row].is_dead()) { |
| continue; |
| } |
| rp = &A[Row[row].start]; |
| rp_end = rp + Row[row].length; |
| while (rp < rp_end) { |
| /* get a column */ |
| col = *rp++; |
| /* add the column, if alive and untagged */ |
| col_thickness = Col[col].shared1.thickness; |
| if (col_thickness > 0 && Col[col].is_alive()) { |
| /* tag column in pivot row */ |
| Col[col].shared1.thickness = -col_thickness; |
| COLAMD_ASSERT(pfree < Alen); |
| /* place column in pivot row */ |
| A[pfree++] = col; |
| pivot_row_degree += col_thickness; |
| } |
| } |
| } |
| |
| /* clear tag on pivot column */ |
| Col[pivot_col].shared1.thickness = pivot_col_thickness; |
| max_deg = numext::maxi(max_deg, pivot_row_degree); |
| |
| /* === Kill all rows used to construct pivot row ==================== */ |
| |
| /* also kill pivot row, temporarily */ |
| cp = &A[Col[pivot_col].start]; |
| cp_end = cp + Col[pivot_col].length; |
| while (cp < cp_end) { |
| /* may be killing an already dead row */ |
| row = *cp++; |
| COLAMD_DEBUG3(("Kill row in pivot col: %d\n", row)); |
| Row[row].kill(); |
| } |
| |
| /* === Select a row index to use as the new pivot row =============== */ |
| |
| pivot_row_length = pfree - pivot_row_start; |
| if (pivot_row_length > 0) { |
| /* pick the "pivot" row arbitrarily (first row in col) */ |
| pivot_row = A[Col[pivot_col].start]; |
| COLAMD_DEBUG3(("Pivotal row is %d\n", pivot_row)); |
| } else { |
| /* there is no pivot row, since it is of zero length */ |
| pivot_row = Empty; |
| COLAMD_ASSERT(pivot_row_length == 0); |
| } |
| COLAMD_ASSERT(Col[pivot_col].length > 0 || pivot_row_length == 0); |
| |
| /* === Approximate degree computation =============================== */ |
| |
| /* Here begins the computation of the approximate degree. The column */ |
| /* score is the sum of the pivot row "length", plus the size of the */ |
| /* set differences of each row in the column minus the pattern of the */ |
| /* pivot row itself. The column ("thickness") itself is also */ |
| /* excluded from the column score (we thus use an approximate */ |
| /* external degree). */ |
| |
| /* The time taken by the following code (compute set differences, and */ |
| /* add them up) is proportional to the size of the data structure */ |
| /* being scanned - that is, the sum of the sizes of each column in */ |
| /* the pivot row. Thus, the amortized time to compute a column score */ |
| /* is proportional to the size of that column (where size, in this */ |
| /* context, is the column "length", or the number of row indices */ |
| /* in that column). The number of row indices in a column is */ |
| /* monotonically non-decreasing, from the length of the original */ |
| /* column on input to colamd. */ |
| |
| /* === Compute set differences ====================================== */ |
| |
| COLAMD_DEBUG3(("** Computing set differences phase. **\n")); |
| |
| /* pivot row is currently dead - it will be revived later. */ |
| |
| COLAMD_DEBUG3(("Pivot row: ")); |
| /* for each column in pivot row */ |
| rp = &A[pivot_row_start]; |
| rp_end = rp + pivot_row_length; |
| while (rp < rp_end) { |
| col = *rp++; |
| COLAMD_ASSERT(Col[col].is_alive() && col != pivot_col); |
| COLAMD_DEBUG3(("Col: %d\n", col)); |
| |
| /* clear tags used to construct pivot row pattern */ |
| col_thickness = -Col[col].shared1.thickness; |
| COLAMD_ASSERT(col_thickness > 0); |
| Col[col].shared1.thickness = col_thickness; |
| |
| /* === Remove column from degree list =========================== */ |
| |
| cur_score = Col[col].shared2.score; |
| prev_col = Col[col].shared3.prev; |
| next_col = Col[col].shared4.degree_next; |
| COLAMD_ASSERT(cur_score >= 0); |
| COLAMD_ASSERT(cur_score <= n_col); |
| COLAMD_ASSERT(cur_score >= Empty); |
| if (prev_col == Empty) { |
| head[cur_score] = next_col; |
| } else { |
| Col[prev_col].shared4.degree_next = next_col; |
| } |
| if (next_col != Empty) { |
| Col[next_col].shared3.prev = prev_col; |
| } |
| |
| /* === Scan the column ========================================== */ |
| |
| cp = &A[Col[col].start]; |
| cp_end = cp + Col[col].length; |
| while (cp < cp_end) { |
| /* get a row */ |
| row = *cp++; |
| /* skip if dead */ |
| if (Row[row].is_dead()) { |
| continue; |
| } |
| row_mark = Row[row].shared2.mark; |
| COLAMD_ASSERT(row != pivot_row); |
| set_difference = row_mark - tag_mark; |
| /* check if the row has been seen yet */ |
| if (set_difference < 0) { |
| COLAMD_ASSERT(Row[row].shared1.degree <= max_deg); |
| set_difference = Row[row].shared1.degree; |
| } |
| /* subtract column thickness from this row's set difference */ |
| set_difference -= col_thickness; |
| COLAMD_ASSERT(set_difference >= 0); |
| /* absorb this row if the set difference becomes zero */ |
| if (set_difference == 0) { |
| COLAMD_DEBUG3(("aggressive absorption. Row: %d\n", row)); |
| Row[row].kill(); |
| } else { |
| /* save the new mark */ |
| Row[row].shared2.mark = set_difference + tag_mark; |
| } |
| } |
| } |
| |
| /* === Add up set differences for each column ======================= */ |
| |
| COLAMD_DEBUG3(("** Adding set differences phase. **\n")); |
| |
| /* for each column in pivot row */ |
| rp = &A[pivot_row_start]; |
| rp_end = rp + pivot_row_length; |
| while (rp < rp_end) { |
| /* get a column */ |
| col = *rp++; |
| COLAMD_ASSERT(Col[col].is_alive() && col != pivot_col); |
| hash = 0; |
| cur_score = 0; |
| cp = &A[Col[col].start]; |
| /* compact the column */ |
| new_cp = cp; |
| cp_end = cp + Col[col].length; |
| |
| COLAMD_DEBUG4(("Adding set diffs for Col: %d.\n", col)); |
| |
| while (cp < cp_end) { |
| /* get a row */ |
| row = *cp++; |
| COLAMD_ASSERT(row >= 0 && row < n_row); |
| /* skip if dead */ |
| if (Row[row].is_dead()) { |
| continue; |
| } |
| row_mark = Row[row].shared2.mark; |
| COLAMD_ASSERT(row_mark > tag_mark); |
| /* compact the column */ |
| *new_cp++ = row; |
| /* compute hash function */ |
| hash += row; |
| /* add set difference */ |
| cur_score += row_mark - tag_mark; |
| /* integer overflow... */ |
| cur_score = numext::mini(cur_score, n_col); |
| } |
| |
| /* recompute the column's length */ |
| Col[col].length = (IndexType)(new_cp - &A[Col[col].start]); |
| |
| /* === Further mass elimination ================================= */ |
| |
| if (Col[col].length == 0) { |
| COLAMD_DEBUG4(("further mass elimination. Col: %d\n", col)); |
| /* nothing left but the pivot row in this column */ |
| Col[col].kill_principal(); |
| pivot_row_degree -= Col[col].shared1.thickness; |
| COLAMD_ASSERT(pivot_row_degree >= 0); |
| /* order it */ |
| Col[col].shared2.order = k; |
| /* increment order count by column thickness */ |
| k += Col[col].shared1.thickness; |
| } else { |
| /* === Prepare for supercolumn detection ==================== */ |
| |
| COLAMD_DEBUG4(("Preparing supercol detection for Col: %d.\n", col)); |
| |
| /* save score so far */ |
| Col[col].shared2.score = cur_score; |
| |
| /* add column to hash table, for supercolumn detection */ |
| hash %= n_col + 1; |
| |
| COLAMD_DEBUG4((" Hash = %d, n_col = %d.\n", hash, n_col)); |
| COLAMD_ASSERT(hash <= n_col); |
| |
| head_column = head[hash]; |
| if (head_column > Empty) { |
| /* degree list "hash" is non-empty, use prev (shared3) of */ |
| /* first column in degree list as head of hash bucket */ |
| first_col = Col[head_column].shared3.headhash; |
| Col[head_column].shared3.headhash = col; |
| } else { |
| /* degree list "hash" is empty, use head as hash bucket */ |
| first_col = -(head_column + 2); |
| head[hash] = -(col + 2); |
| } |
| Col[col].shared4.hash_next = first_col; |
| |
| /* save hash function in Col [col].shared3.hash */ |
| Col[col].shared3.hash = (IndexType)hash; |
| COLAMD_ASSERT(Col[col].is_alive()); |
| } |
| } |
| |
| /* The approximate external column degree is now computed. */ |
| |
| /* === Supercolumn detection ======================================== */ |
| |
| COLAMD_DEBUG3(("** Supercolumn detection phase. **\n")); |
| |
| Colamd::detect_super_cols(Col, A, head, pivot_row_start, pivot_row_length); |
| |
| /* === Kill the pivotal column ====================================== */ |
| |
| Col[pivot_col].kill_principal(); |
| |
| /* === Clear mark =================================================== */ |
| |
| tag_mark += (max_deg + 1); |
| if (tag_mark >= max_mark) { |
| COLAMD_DEBUG2(("clearing tag_mark\n")); |
| tag_mark = Colamd::clear_mark(n_row, Row); |
| } |
| |
| /* === Finalize the new pivot row, and column scores ================ */ |
| |
| COLAMD_DEBUG3(("** Finalize scores phase. **\n")); |
| |
| /* for each column in pivot row */ |
| rp = &A[pivot_row_start]; |
| /* compact the pivot row */ |
| new_rp = rp; |
| rp_end = rp + pivot_row_length; |
| while (rp < rp_end) { |
| col = *rp++; |
| /* skip dead columns */ |
| if (Col[col].is_dead()) { |
| continue; |
| } |
| *new_rp++ = col; |
| /* add new pivot row to column */ |
| A[Col[col].start + (Col[col].length++)] = pivot_row; |
| |
| /* retrieve score so far and add on pivot row's degree. */ |
| /* (we wait until here for this in case the pivot */ |
| /* row's degree was reduced due to mass elimination). */ |
| cur_score = Col[col].shared2.score + pivot_row_degree; |
| |
| /* calculate the max possible score as the number of */ |
| /* external columns minus the 'k' value minus the */ |
| /* columns thickness */ |
| max_score = n_col - k - Col[col].shared1.thickness; |
| |
| /* make the score the external degree of the union-of-rows */ |
| cur_score -= Col[col].shared1.thickness; |
| |
| /* make sure score is less or equal than the max score */ |
| cur_score = numext::mini(cur_score, max_score); |
| COLAMD_ASSERT(cur_score >= 0); |
| |
| /* store updated score */ |
| Col[col].shared2.score = cur_score; |
| |
| /* === Place column back in degree list ========================= */ |
| |
| COLAMD_ASSERT(min_score >= 0); |
| COLAMD_ASSERT(min_score <= n_col); |
| COLAMD_ASSERT(cur_score >= 0); |
| COLAMD_ASSERT(cur_score <= n_col); |
| COLAMD_ASSERT(head[cur_score] >= Empty); |
| next_col = head[cur_score]; |
| Col[col].shared4.degree_next = next_col; |
| Col[col].shared3.prev = Empty; |
| if (next_col != Empty) { |
| Col[next_col].shared3.prev = col; |
| } |
| head[cur_score] = col; |
| |
| /* see if this score is less than current min */ |
| min_score = numext::mini(min_score, cur_score); |
| } |
| |
| /* === Resurrect the new pivot row ================================== */ |
| |
| if (pivot_row_degree > 0) { |
| /* update pivot row length to reflect any cols that were killed */ |
| /* during super-col detection and mass elimination */ |
| Row[pivot_row].start = pivot_row_start; |
| Row[pivot_row].length = (IndexType)(new_rp - &A[pivot_row_start]); |
| Row[pivot_row].shared1.degree = pivot_row_degree; |
| Row[pivot_row].shared2.mark = 0; |
| /* pivot row is no longer dead */ |
| } |
| } |
| |
| /* === All principal columns have now been ordered ====================== */ |
| |
| return (ngarbage); |
| } |
| |
| /* ========================================================================== */ |
| /* === order_children ======================================================= */ |
| /* ========================================================================== */ |
| |
| /* |
| The find_ordering routine has ordered all of the principal columns (the |
| representatives of the supercolumns). The non-principal columns have not |
| yet been ordered. This routine orders those columns by walking up the |
| parent tree (a column is a child of the column which absorbed it). The |
| final permutation vector is then placed in p [0 ... n_col-1], with p [0] |
| being the first column, and p [n_col-1] being the last. It doesn't look |
| like it at first glance, but be assured that this routine takes time linear |
| in the number of columns. Although not immediately obvious, the time |
| taken by this routine is O (n_col), that is, linear in the number of |
| columns. Not user-callable. |
| */ |
| template <typename IndexType> |
| static inline void order_children( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_col, /* number of columns of A */ |
| ColStructure<IndexType> Col[], /* of size n_col+1 */ |
| IndexType p[] /* p [0 ... n_col-1] is the column permutation*/ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType i; /* loop counter for all columns */ |
| IndexType c; /* column index */ |
| IndexType parent; /* index of column's parent */ |
| IndexType order; /* column's order */ |
| |
| /* === Order each non-principal column ================================== */ |
| |
| for (i = 0; i < n_col; i++) { |
| /* find an un-ordered non-principal column */ |
| COLAMD_ASSERT(col_is_dead(Col, i)); |
| if (!Col[i].is_dead_principal() && Col[i].shared2.order == Empty) { |
| parent = i; |
| /* once found, find its principal parent */ |
| do { |
| parent = Col[parent].shared1.parent; |
| } while (!Col[parent].is_dead_principal()); |
| |
| /* now, order all un-ordered non-principal columns along path */ |
| /* to this parent. collapse tree at the same time */ |
| c = i; |
| /* get order of parent */ |
| order = Col[parent].shared2.order; |
| |
| do { |
| COLAMD_ASSERT(Col[c].shared2.order == Empty); |
| |
| /* order this column */ |
| Col[c].shared2.order = order++; |
| /* collapse tree */ |
| Col[c].shared1.parent = parent; |
| |
| /* get immediate parent of this column */ |
| c = Col[c].shared1.parent; |
| |
| /* continue until we hit an ordered column. There are */ |
| /* guaranteed not to be anymore unordered columns */ |
| /* above an ordered column */ |
| } while (Col[c].shared2.order == Empty); |
| |
| /* re-order the super_col parent to largest order for this group */ |
| Col[parent].shared2.order = order; |
| } |
| } |
| |
| /* === Generate the permutation ========================================= */ |
| |
| for (c = 0; c < n_col; c++) { |
| p[Col[c].shared2.order] = c; |
| } |
| } |
| |
| /* ========================================================================== */ |
| /* === detect_super_cols ==================================================== */ |
| /* ========================================================================== */ |
| |
| /* |
| Detects supercolumns by finding matches between columns in the hash buckets. |
| Check amongst columns in the set A [row_start ... row_start + row_length-1]. |
| The columns under consideration are currently *not* in the degree lists, |
| and have already been placed in the hash buckets. |
| |
| The hash bucket for columns whose hash function is equal to h is stored |
| as follows: |
| |
| if head [h] is >= 0, then head [h] contains a degree list, so: |
| |
| head [h] is the first column in degree bucket h. |
| Col [head [h]].headhash gives the first column in hash bucket h. |
| |
| otherwise, the degree list is empty, and: |
| |
| -(head [h] + 2) is the first column in hash bucket h. |
| |
| For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous |
| column" pointer. Col [c].shared3.hash is used instead as the hash number |
| for that column. The value of Col [c].shared4.hash_next is the next column |
| in the same hash bucket. |
| |
| Assuming no, or "few" hash collisions, the time taken by this routine is |
| linear in the sum of the sizes (lengths) of each column whose score has |
| just been computed in the approximate degree computation. |
| Not user-callable. |
| */ |
| template <typename IndexType> |
| static void detect_super_cols( |
| /* === Parameters ======================================================= */ |
| |
| ColStructure<IndexType> Col[], /* of size n_col+1 */ |
| IndexType A[], /* row indices of A */ |
| IndexType head[], /* head of degree lists and hash buckets */ |
| IndexType row_start, /* pointer to set of columns to check */ |
| IndexType row_length /* number of columns to check */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType hash; /* hash value for a column */ |
| IndexType *rp; /* pointer to a row */ |
| IndexType c; /* a column index */ |
| IndexType super_c; /* column index of the column to absorb into */ |
| IndexType *cp1; /* column pointer for column super_c */ |
| IndexType *cp2; /* column pointer for column c */ |
| IndexType length; /* length of column super_c */ |
| IndexType prev_c; /* column preceding c in hash bucket */ |
| IndexType i; /* loop counter */ |
| IndexType *rp_end; /* pointer to the end of the row */ |
| IndexType col; /* a column index in the row to check */ |
| IndexType head_column; /* first column in hash bucket or degree list */ |
| IndexType first_col; /* first column in hash bucket */ |
| |
| /* === Consider each column in the row ================================== */ |
| |
| rp = &A[row_start]; |
| rp_end = rp + row_length; |
| while (rp < rp_end) { |
| col = *rp++; |
| if (Col[col].is_dead()) { |
| continue; |
| } |
| |
| /* get hash number for this column */ |
| hash = Col[col].shared3.hash; |
| COLAMD_ASSERT(hash <= n_col); |
| |
| /* === Get the first column in this hash bucket ===================== */ |
| |
| head_column = head[hash]; |
| if (head_column > Empty) { |
| first_col = Col[head_column].shared3.headhash; |
| } else { |
| first_col = -(head_column + 2); |
| } |
| |
| /* === Consider each column in the hash bucket ====================== */ |
| |
| for (super_c = first_col; super_c != Empty; super_c = Col[super_c].shared4.hash_next) { |
| COLAMD_ASSERT(Col[super_c].is_alive()); |
| COLAMD_ASSERT(Col[super_c].shared3.hash == hash); |
| length = Col[super_c].length; |
| |
| /* prev_c is the column preceding column c in the hash bucket */ |
| prev_c = super_c; |
| |
| /* === Compare super_c with all columns after it ================ */ |
| |
| for (c = Col[super_c].shared4.hash_next; c != Empty; c = Col[c].shared4.hash_next) { |
| COLAMD_ASSERT(c != super_c); |
| COLAMD_ASSERT(Col[c].is_alive()); |
| COLAMD_ASSERT(Col[c].shared3.hash == hash); |
| |
| /* not identical if lengths or scores are different */ |
| if (Col[c].length != length || Col[c].shared2.score != Col[super_c].shared2.score) { |
| prev_c = c; |
| continue; |
| } |
| |
| /* compare the two columns */ |
| cp1 = &A[Col[super_c].start]; |
| cp2 = &A[Col[c].start]; |
| |
| for (i = 0; i < length; i++) { |
| /* the columns are "clean" (no dead rows) */ |
| COLAMD_ASSERT(cp1->is_alive()); |
| COLAMD_ASSERT(cp2->is_alive()); |
| /* row indices will same order for both supercols, */ |
| /* no gather scatter necessary */ |
| if (*cp1++ != *cp2++) { |
| break; |
| } |
| } |
| |
| /* the two columns are different if the for-loop "broke" */ |
| if (i != length) { |
| prev_c = c; |
| continue; |
| } |
| |
| /* === Got it! two columns are identical =================== */ |
| |
| COLAMD_ASSERT(Col[c].shared2.score == Col[super_c].shared2.score); |
| |
| Col[super_c].shared1.thickness += Col[c].shared1.thickness; |
| Col[c].shared1.parent = super_c; |
| Col[c].kill_non_principal(); |
| /* order c later, in order_children() */ |
| Col[c].shared2.order = Empty; |
| /* remove c from hash bucket */ |
| Col[prev_c].shared4.hash_next = Col[c].shared4.hash_next; |
| } |
| } |
| |
| /* === Empty this hash bucket ======================================= */ |
| |
| if (head_column > Empty) { |
| /* corresponding degree list "hash" is not empty */ |
| Col[head_column].shared3.headhash = Empty; |
| } else { |
| /* corresponding degree list "hash" is empty */ |
| head[hash] = Empty; |
| } |
| } |
| } |
| |
| /* ========================================================================== */ |
| /* === garbage_collection =================================================== */ |
| /* ========================================================================== */ |
| |
| /* |
| Defragments and compacts columns and rows in the workspace A. Used when |
| all available memory has been used while performing row merging. Returns |
| the index of the first free position in A, after garbage collection. The |
| time taken by this routine is linear is the size of the array A, which is |
| itself linear in the number of nonzeros in the input matrix. |
| Not user-callable. |
| */ |
| template <typename IndexType> |
| static IndexType garbage_collection /* returns the new value of pfree */ |
| ( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_row, /* number of rows */ |
| IndexType n_col, /* number of columns */ |
| RowStructure<IndexType> Row[], /* row info */ |
| ColStructure<IndexType> Col[], /* column info */ |
| IndexType A[], /* A [0 ... Alen-1] holds the matrix */ |
| IndexType *pfree /* &A [0] ... pfree is in use */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType *psrc; /* source pointer */ |
| IndexType *pdest; /* destination pointer */ |
| IndexType j; /* counter */ |
| IndexType r; /* a row index */ |
| IndexType c; /* a column index */ |
| IndexType length; /* length of a row or column */ |
| |
| /* === Defragment the columns =========================================== */ |
| |
| pdest = &A[0]; |
| for (c = 0; c < n_col; c++) { |
| if (Col[c].is_alive()) { |
| psrc = &A[Col[c].start]; |
| |
| /* move and compact the column */ |
| COLAMD_ASSERT(pdest <= psrc); |
| Col[c].start = (IndexType)(pdest - &A[0]); |
| length = Col[c].length; |
| for (j = 0; j < length; j++) { |
| r = *psrc++; |
| if (Row[r].is_alive()) { |
| *pdest++ = r; |
| } |
| } |
| Col[c].length = (IndexType)(pdest - &A[Col[c].start]); |
| } |
| } |
| |
| /* === Prepare to defragment the rows =================================== */ |
| |
| for (r = 0; r < n_row; r++) { |
| if (Row[r].is_alive()) { |
| if (Row[r].length == 0) { |
| /* this row is of zero length. cannot compact it, so kill it */ |
| COLAMD_DEBUG3(("Defrag row kill\n")); |
| Row[r].kill(); |
| } else { |
| /* save first column index in Row [r].shared2.first_column */ |
| psrc = &A[Row[r].start]; |
| Row[r].shared2.first_column = *psrc; |
| COLAMD_ASSERT(Row[r].is_alive()); |
| /* flag the start of the row with the one's complement of row */ |
| *psrc = ones_complement(r); |
| } |
| } |
| } |
| |
| /* === Defragment the rows ============================================== */ |
| |
| psrc = pdest; |
| while (psrc < pfree) { |
| /* find a negative number ... the start of a row */ |
| if (*psrc++ < 0) { |
| psrc--; |
| /* get the row index */ |
| r = ones_complement(*psrc); |
| COLAMD_ASSERT(r >= 0 && r < n_row); |
| /* restore first column index */ |
| *psrc = Row[r].shared2.first_column; |
| COLAMD_ASSERT(Row[r].is_alive()); |
| |
| /* move and compact the row */ |
| COLAMD_ASSERT(pdest <= psrc); |
| Row[r].start = (IndexType)(pdest - &A[0]); |
| length = Row[r].length; |
| for (j = 0; j < length; j++) { |
| c = *psrc++; |
| if (Col[c].is_alive()) { |
| *pdest++ = c; |
| } |
| } |
| Row[r].length = (IndexType)(pdest - &A[Row[r].start]); |
| } |
| } |
| /* ensure we found all the rows */ |
| COLAMD_ASSERT(debug_rows == 0); |
| |
| /* === Return the new value of pfree ==================================== */ |
| |
| return ((IndexType)(pdest - &A[0])); |
| } |
| |
| /* ========================================================================== */ |
| /* === clear_mark =========================================================== */ |
| /* ========================================================================== */ |
| |
| /* |
| Clears the Row [].shared2.mark array, and returns the new tag_mark. |
| Return value is the new tag_mark. Not user-callable. |
| */ |
| template <typename IndexType> |
| static inline IndexType clear_mark /* return the new value for tag_mark */ |
| ( |
| /* === Parameters ======================================================= */ |
| |
| IndexType n_row, /* number of rows in A */ |
| RowStructure<IndexType> Row[] /* Row [0 ... n_row-1].shared2.mark is set to zero */ |
| ) { |
| /* === Local variables ================================================== */ |
| |
| IndexType r; |
| |
| for (r = 0; r < n_row; r++) { |
| if (Row[r].is_alive()) { |
| Row[r].shared2.mark = 0; |
| } |
| } |
| return (1); |
| } |
| |
| } // namespace Colamd |
| |
| } // namespace internal |
| #endif |
