diff options
Diffstat (limited to 'eigen/Eigen/src/OrderingMethods')
| -rw-r--r-- | eigen/Eigen/src/OrderingMethods/Amd.h | 83 | ||||
| -rw-r--r-- | eigen/Eigen/src/OrderingMethods/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h | 412 | ||||
| -rw-r--r-- | eigen/Eigen/src/OrderingMethods/Ordering.h | 51 |
4 files changed, 275 insertions, 277 deletions
diff --git a/eigen/Eigen/src/OrderingMethods/Amd.h b/eigen/Eigen/src/OrderingMethods/Amd.h index 658b954..f91ecb2 100644 --- a/eigen/Eigen/src/OrderingMethods/Amd.h +++ b/eigen/Eigen/src/OrderingMethods/Amd.h @@ -41,10 +41,10 @@ template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); } /* clear w */ -template<typename Index> -static int cs_wclear (Index mark, Index lemax, Index *w, Index n) +template<typename StorageIndex> +static StorageIndex cs_wclear (StorageIndex mark, StorageIndex lemax, StorageIndex *w, StorageIndex n) { - Index k; + StorageIndex k; if(mark < 2 || (mark + lemax < 0)) { for(k = 0; k < n; k++) @@ -56,10 +56,10 @@ static int cs_wclear (Index mark, Index lemax, Index *w, Index n) } /* depth-first search and postorder of a tree rooted at node j */ -template<typename Index> -Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack) +template<typename StorageIndex> +StorageIndex cs_tdfs(StorageIndex j, StorageIndex k, StorageIndex *head, const StorageIndex *next, StorageIndex *post, StorageIndex *stack) { - int i, p, top = 0; + StorageIndex i, p, top = 0; if(!head || !next || !post || !stack) return (-1); /* check inputs */ stack[0] = j; /* place j on the stack */ while (top >= 0) /* while (stack is not empty) */ @@ -84,42 +84,45 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind /** \internal * \ingroup OrderingMethods_Module * Approximate minimum degree ordering algorithm. - * \returns the permutation P reducing the fill-in of the input matrix \a C - * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional. + * + * \param[in] C the input selfadjoint matrix stored in compressed column major format. + * \param[out] perm the permutation P reducing the fill-in of the input matrix \a C + * + * Note that the input matrix \a C must be complete, that is both the upper and lower parts have to be stored, as well as the diagonal entries. * On exit the values of C are destroyed */ -template<typename Scalar, typename Index> -void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm) +template<typename Scalar, typename StorageIndex> +void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,StorageIndex>& C, PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm) { using std::sqrt; - int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, - k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, - ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t; - unsigned int h; + StorageIndex d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, + k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, + ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t, h; - Index n = C.cols(); - dense = std::max<Index> (16, Index(10 * sqrt(double(n)))); /* find dense threshold */ - dense = std::min<Index> (n-2, dense); + StorageIndex n = StorageIndex(C.cols()); + dense = std::max<StorageIndex> (16, StorageIndex(10 * sqrt(double(n)))); /* find dense threshold */ + dense = (std::min)(n-2, dense); - Index cnz = C.nonZeros(); + StorageIndex cnz = StorageIndex(C.nonZeros()); perm.resize(n+1); t = cnz + cnz/5 + 2*n; /* add elbow room to C */ C.resizeNonZeros(t); - Index* W = new Index[8*(n+1)]; /* get workspace */ - Index* len = W; - Index* nv = W + (n+1); - Index* next = W + 2*(n+1); - Index* head = W + 3*(n+1); - Index* elen = W + 4*(n+1); - Index* degree = W + 5*(n+1); - Index* w = W + 6*(n+1); - Index* hhead = W + 7*(n+1); - Index* last = perm.indices().data(); /* use P as workspace for last */ + // get workspace + ei_declare_aligned_stack_constructed_variable(StorageIndex,W,8*(n+1),0); + StorageIndex* len = W; + StorageIndex* nv = W + (n+1); + StorageIndex* next = W + 2*(n+1); + StorageIndex* head = W + 3*(n+1); + StorageIndex* elen = W + 4*(n+1); + StorageIndex* degree = W + 5*(n+1); + StorageIndex* w = W + 6*(n+1); + StorageIndex* hhead = W + 7*(n+1); + StorageIndex* last = perm.indices().data(); /* use P as workspace for last */ /* --- Initialize quotient graph ---------------------------------------- */ - Index* Cp = C.outerIndexPtr(); - Index* Ci = C.innerIndexPtr(); + StorageIndex* Cp = C.outerIndexPtr(); + StorageIndex* Ci = C.innerIndexPtr(); for(k = 0; k < n; k++) len[k] = Cp[k+1] - Cp[k]; len[n] = 0; @@ -136,7 +139,7 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation elen[i] = 0; // Ek of node i is empty degree[i] = len[i]; // degree of node i } - mark = internal::cs_wclear<Index>(0, 0, w, n); /* clear w */ + mark = internal::cs_wclear<StorageIndex>(0, 0, w, n); /* clear w */ /* --- Initialize degree lists ------------------------------------------ */ for(i = 0; i < n; i++) @@ -260,7 +263,7 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation elen[k] = -2; /* k is now an element */ /* --- Find set differences ----------------------------------------- */ - mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */ + mark = internal::cs_wclear<StorageIndex>(mark, lemax, w, n); /* clear w if necessary */ for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */ { i = Ci[pk]; @@ -330,7 +333,7 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation } else { - degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */ + degree[i] = std::min<StorageIndex> (degree[i], d); /* update degree(i) */ Ci[pn] = Ci[p3]; /* move first node to end */ Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */ Ci[p1] = k; /* add k as 1st element in of Ei */ @@ -338,12 +341,12 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation h %= n; /* finalize hash of i */ next[i] = hhead[h]; /* place i in hash bucket */ hhead[h] = i; - last[i] = h; /* save hash of i in last[i] */ + last[i] = h; /* save hash of i in last[i] */ } } /* scan2 is done */ degree[k] = dk; /* finalize |Lk| */ - lemax = std::max<Index>(lemax, dk); - mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */ + lemax = std::max<StorageIndex>(lemax, dk); + mark = internal::cs_wclear<StorageIndex>(mark+lemax, lemax, w, n); /* clear w */ /* --- Supernode detection ------------------------------------------ */ for(pk = pk1; pk < pk2; pk++) @@ -391,12 +394,12 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */ nv[i] = nvi; /* restore nv[i] */ d = degree[i] + dk - nvi; /* compute external degree(i) */ - d = std::min<Index> (d, n - nel - nvi); + d = std::min<StorageIndex> (d, n - nel - nvi); if(head[d] != -1) last[head[d]] = i; next[i] = head[d]; /* put i back in degree list */ last[i] = -1; head[d] = i; - mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */ + mindeg = std::min<StorageIndex> (mindeg, d); /* find new minimum degree */ degree[i] = d; Ci[p++] = i; /* place i in Lk */ } @@ -429,12 +432,10 @@ void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, Permutation } for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */ { - if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w); + if(Cp[i] == -1) k = internal::cs_tdfs<StorageIndex>(i, k, head, next, perm.indices().data(), w); } perm.indices().conservativeResize(n); - - delete[] W; } } // namespace internal diff --git a/eigen/Eigen/src/OrderingMethods/CMakeLists.txt b/eigen/Eigen/src/OrderingMethods/CMakeLists.txt deleted file mode 100644 index 9f4bb27..0000000 --- a/eigen/Eigen/src/OrderingMethods/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_OrderingMethods_SRCS "*.h") - -INSTALL(FILES - ${Eigen_OrderingMethods_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/OrderingMethods COMPONENT Devel - ) diff --git a/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h b/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h index 359fd44..da85b4d 100644 --- a/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h +++ b/eigen/Eigen/src/OrderingMethods/Eigen_Colamd.h @@ -128,54 +128,54 @@ namespace internal { /* ========================================================================== */ // == Row and Column structures == -template <typename Index> +template <typename IndexType> struct colamd_col { - Index start ; /* index for A of first row in this column, or DEAD */ + IndexType start ; /* index for A of first row in this column, or DEAD */ /* if column is dead */ - Index length ; /* number of rows in this column */ + IndexType length ; /* number of rows in this column */ union { - Index thickness ; /* number of original columns represented by this */ + IndexType thickness ; /* number of original columns represented by this */ /* col, if the column is alive */ - Index parent ; /* parent in parent tree super-column structure, if */ + IndexType parent ; /* parent in parent tree super-column structure, if */ /* the column is dead */ } shared1 ; union { - Index score ; /* the score used to maintain heap, if col is alive */ - Index order ; /* pivot ordering of this column, if col is dead */ + 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 { - Index headhash ; /* head of a hash bucket, if col is at the head of */ + IndexType headhash ; /* head of a hash bucket, if col is at the head of */ /* a degree list */ - Index hash ; /* hash value, if col is not in a degree list */ - Index prev ; /* previous column in degree list, if col is in a */ + 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 { - Index degree_next ; /* next column, if col is in a degree list */ - Index hash_next ; /* next column, if col is in a hash list */ + 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 ; }; -template <typename Index> +template <typename IndexType> struct Colamd_Row { - Index start ; /* index for A of first col in this row */ - Index length ; /* number of principal columns in this row */ + IndexType start ; /* index for A of first col in this row */ + IndexType length ; /* number of principal columns in this row */ union { - Index degree ; /* number of principal & non-principal columns in row */ - Index p ; /* used as a row pointer in init_rows_cols () */ + IndexType degree ; /* number of principal & non-principal columns in row */ + IndexType p ; /* used as a row pointer in init_rows_cols () */ } shared1 ; union { - Index mark ; /* for computing set differences and marking dead rows*/ - Index first_column ;/* first column in row (used in garbage collection) */ + IndexType mark ; /* for computing set differences and marking dead rows*/ + IndexType first_column ;/* first column in row (used in garbage collection) */ } shared2 ; }; @@ -195,38 +195,38 @@ struct Colamd_Row This macro is not needed when using symamd. - Explicit typecast to Index added Sept. 23, 2002, COLAMD version 2.2, to avoid + Explicit typecast to IndexType added Sept. 23, 2002, COLAMD version 2.2, to avoid gcc -pedantic warning messages. */ -template <typename Index> -inline Index colamd_c(Index n_col) -{ return Index( ((n_col) + 1) * sizeof (colamd_col<Index>) / sizeof (Index) ) ; } +template <typename IndexType> +inline IndexType colamd_c(IndexType n_col) +{ return IndexType( ((n_col) + 1) * sizeof (colamd_col<IndexType>) / sizeof (IndexType) ) ; } -template <typename Index> -inline Index colamd_r(Index n_row) -{ return Index(((n_row) + 1) * sizeof (Colamd_Row<Index>) / sizeof (Index)); } +template <typename IndexType> +inline IndexType colamd_r(IndexType n_row) +{ return IndexType(((n_row) + 1) * sizeof (Colamd_Row<IndexType>) / sizeof (IndexType)); } // Prototypes of non-user callable routines -template <typename Index> -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] ); +template <typename IndexType> +static IndexType init_rows_cols (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> col [], IndexType A [], IndexType p [], IndexType stats[COLAMD_STATS] ); -template <typename Index> -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); +template <typename IndexType> +static void init_scoring (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], double knobs[COLAMD_KNOBS], IndexType *p_n_row2, IndexType *p_n_col2, IndexType *p_max_deg); -template <typename Index> -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); +template <typename IndexType> +static IndexType find_ordering (IndexType n_row, IndexType n_col, IndexType Alen, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType n_col2, IndexType max_deg, IndexType pfree); -template <typename Index> -static void order_children (Index n_col, colamd_col<Index> Col [], Index p []); +template <typename IndexType> +static void order_children (IndexType n_col, colamd_col<IndexType> Col [], IndexType p []); -template <typename Index> -static void detect_super_cols (colamd_col<Index> Col [], Index A [], Index head [], Index row_start, Index row_length ) ; +template <typename IndexType> +static void detect_super_cols (colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType row_start, IndexType row_length ) ; -template <typename Index> -static Index garbage_collection (Index n_row, Index n_col, Colamd_Row<Index> Row [], colamd_col<Index> Col [], Index A [], Index *pfree) ; +template <typename IndexType> +static IndexType garbage_collection (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType *pfree) ; -template <typename Index> -static inline Index clear_mark (Index n_row, Colamd_Row<Index> Row [] ) ; +template <typename IndexType> +static inline IndexType clear_mark (IndexType n_row, Colamd_Row<IndexType> Row [] ) ; /* === No debugging ========================================================= */ @@ -253,8 +253,8 @@ static inline Index clear_mark (Index n_row, Colamd_Row<Index> Row [] ) ; * \param n_col number of columns in A * \return recommended value of Alen for use by colamd */ -template <typename Index> -inline Index colamd_recommended ( Index nnz, Index n_row, Index n_col) +template <typename IndexType> +inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col) { if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) return (-1); @@ -318,22 +318,22 @@ static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS]) * \param knobs parameter settings for colamd * \param stats colamd output statistics and error codes */ -template <typename Index> -static bool colamd(Index n_row, Index n_col, Index Alen, Index *A, Index *p, double knobs[COLAMD_KNOBS], Index stats[COLAMD_STATS]) +template <typename IndexType> +static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p, double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS]) { /* === Local variables ================================================== */ - Index i ; /* loop index */ - Index nnz ; /* nonzeros in A */ - Index Row_size ; /* size of Row [], in integers */ - Index Col_size ; /* size of Col [], in integers */ - Index need ; /* minimum required length of A */ - Colamd_Row<Index> *Row ; /* pointer into A of Row [0..n_row] array */ - colamd_col<Index> *Col ; /* pointer into A of Col [0..n_col] array */ - Index n_col2 ; /* number of non-dense, non-empty columns */ - Index n_row2 ; /* number of non-dense, non-empty rows */ - Index ngarbage ; /* number of garbage collections performed */ - Index max_deg ; /* maximum row degree */ + 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_Row<IndexType> *Row ; /* pointer into A of Row [0..n_row] array */ + colamd_col<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 [COLAMD_KNOBS] ; /* default knobs array */ @@ -424,8 +424,8 @@ static bool colamd(Index n_row, Index n_col, Index Alen, Index *A, Index *p, dou } Alen -= Col_size + Row_size ; - Col = (colamd_col<Index> *) &A [Alen] ; - Row = (Colamd_Row<Index> *) &A [Alen + Col_size] ; + Col = (colamd_col<IndexType> *) &A [Alen] ; + Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ; /* === Construct the row and column data structures ===================== */ @@ -478,29 +478,29 @@ static bool colamd(Index n_row, Index n_col, Index Alen, Index *A, Index *p, dou column form of the matrix. Returns false if the matrix is invalid, true otherwise. Not user-callable. */ -template <typename Index> -static Index init_rows_cols /* returns true if OK, or false otherwise */ +template <typename IndexType> +static IndexType init_rows_cols /* returns true if OK, or false otherwise */ ( /* === Parameters ======================================================= */ - Index n_row, /* number of rows of A */ - Index n_col, /* number of columns of A */ - Colamd_Row<Index> Row [], /* of size n_row+1 */ - colamd_col<Index> Col [], /* of size n_col+1 */ - Index A [], /* row indices of A, of size Alen */ - Index p [], /* pointers to columns in A, of size n_col+1 */ - Index stats [COLAMD_STATS] /* colamd statistics */ + IndexType n_row, /* number of rows of A */ + IndexType n_col, /* number of columns of A */ + Colamd_Row<IndexType> Row [], /* of size n_row+1 */ + colamd_col<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 [COLAMD_STATS] /* colamd statistics */ ) { /* === Local variables ================================================== */ - Index col ; /* a column index */ - Index row ; /* a row index */ - Index *cp ; /* a column pointer */ - Index *cp_end ; /* a pointer to the end of a column */ - Index *rp ; /* a row pointer */ - Index *rp_end ; /* a pointer to the end of a row */ - Index last_row ; /* previous row */ + 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 ==================== */ @@ -694,46 +694,46 @@ static Index init_rows_cols /* returns true if OK, or false otherwise */ 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 Index> +template <typename IndexType> static void init_scoring ( /* === Parameters ======================================================= */ - Index n_row, /* number of rows of A */ - Index n_col, /* number of columns of A */ - Colamd_Row<Index> Row [], /* of size n_row+1 */ - colamd_col<Index> Col [], /* of size n_col+1 */ - Index A [], /* column form and row form of A */ - Index head [], /* of size n_col+1 */ + IndexType n_row, /* number of rows of A */ + IndexType n_col, /* number of columns of A */ + Colamd_Row<IndexType> Row [], /* of size n_row+1 */ + colamd_col<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 [COLAMD_KNOBS],/* parameters */ - Index *p_n_row2, /* number of non-dense, non-empty rows */ - Index *p_n_col2, /* number of non-dense, non-empty columns */ - Index *p_max_deg /* maximum row degree */ + 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 ================================================== */ - Index c ; /* a column index */ - Index r, row ; /* a row index */ - Index *cp ; /* a column pointer */ - Index deg ; /* degree of a row or column */ - Index *cp_end ; /* a pointer to the end of a column */ - Index *new_cp ; /* new column pointer */ - Index col_length ; /* length of pruned column */ - Index score ; /* current column score */ - Index n_col2 ; /* number of non-dense, non-empty columns */ - Index n_row2 ; /* number of non-dense, non-empty rows */ - Index dense_row_count ; /* remove rows with more entries than this */ - Index dense_col_count ; /* remove cols with more entries than this */ - Index min_score ; /* smallest column score */ - Index max_deg ; /* maximum row degree */ - Index next_col ; /* Used to add to degree list.*/ + 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 = std::max<Index>(0, (std::min)(Index(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ; - dense_col_count = std::max<Index>(0, (std::min)(Index(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ; + dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ; + dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ; COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; max_deg = 0 ; n_col2 = n_col ; @@ -797,7 +797,7 @@ static void init_scoring else { /* keep track of max degree of remaining rows */ - max_deg = (std::max)(max_deg, deg) ; + max_deg = numext::maxi(max_deg, deg) ; } } COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; @@ -835,10 +835,10 @@ static void init_scoring /* add row's external degree */ score += Row [row].shared1.degree - 1 ; /* guard against integer overflow */ - score = (std::min)(score, n_col) ; + score = numext::mini(score, n_col) ; } /* determine pruned column length */ - col_length = (Index) (new_cp - &A [Col [c].start]) ; + 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" */ @@ -907,7 +907,7 @@ static void init_scoring head [score] = c ; /* see if this score is less than current min */ - min_score = (std::min)(min_score, score) ; + min_score = numext::mini(min_score, score) ; } @@ -931,56 +931,56 @@ static void init_scoring (no supercolumns on input). Uses a minimum approximate column minimum degree ordering method. Not user-callable. */ -template <typename Index> -static Index find_ordering /* return the number of garbage collections */ +template <typename IndexType> +static IndexType find_ordering /* return the number of garbage collections */ ( /* === Parameters ======================================================= */ - Index n_row, /* number of rows of A */ - Index n_col, /* number of columns of A */ - Index Alen, /* size of A, 2*nnz + n_col or larger */ - Colamd_Row<Index> Row [], /* of size n_row+1 */ - colamd_col<Index> Col [], /* of size n_col+1 */ - Index A [], /* column form and row form of A */ - Index head [], /* of size n_col+1 */ - Index n_col2, /* Remaining columns to order */ - Index max_deg, /* Maximum row degree */ - Index pfree /* index of first free slot (2*nnz on entry) */ + 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 */ + Colamd_Row<IndexType> Row [], /* of size n_row+1 */ + colamd_col<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 ================================================== */ - Index k ; /* current pivot ordering step */ - Index pivot_col ; /* current pivot column */ - Index *cp ; /* a column pointer */ - Index *rp ; /* a row pointer */ - Index pivot_row ; /* current pivot row */ - Index *new_cp ; /* modified column pointer */ - Index *new_rp ; /* modified row pointer */ - Index pivot_row_start ; /* pointer to start of pivot row */ - Index pivot_row_degree ; /* number of columns in pivot row */ - Index pivot_row_length ; /* number of supercolumns in pivot row */ - Index pivot_col_score ; /* score of pivot column */ - Index needed_memory ; /* free space needed for pivot row */ - Index *cp_end ; /* pointer to the end of a column */ - Index *rp_end ; /* pointer to the end of a row */ - Index row ; /* a row index */ - Index col ; /* a column index */ - Index max_score ; /* maximum possible score */ - Index cur_score ; /* score of current column */ + 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 */ - Index head_column ; /* head of hash bucket */ - Index first_col ; /* first column in hash bucket */ - Index tag_mark ; /* marker value for mark array */ - Index row_mark ; /* Row [row].shared2.mark */ - Index set_difference ; /* set difference size of row with pivot row */ - Index min_score ; /* smallest column score */ - Index col_thickness ; /* "thickness" (no. of columns in a supercol) */ - Index max_mark ; /* maximum value of tag_mark */ - Index pivot_col_thickness ; /* number of columns represented by pivot col */ - Index prev_col ; /* Used by Dlist operations. */ - Index next_col ; /* Used by Dlist operations. */ - Index ngarbage ; /* number of garbage collections performed */ + 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 ==================================== */ @@ -1004,7 +1004,7 @@ static Index find_ordering /* return the number of garbage collections */ COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ; /* get pivot column from head of minimum degree list */ - while (head [min_score] == COLAMD_EMPTY && min_score < n_col) + while (min_score < n_col && head [min_score] == COLAMD_EMPTY) { min_score++ ; } @@ -1033,7 +1033,7 @@ static Index find_ordering /* return the number of garbage collections */ /* === Garbage_collection, if necessary ============================= */ - needed_memory = (std::min)(pivot_col_score, n_col - k) ; + needed_memory = numext::mini(pivot_col_score, n_col - k) ; if (pfree + needed_memory >= Alen) { pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; @@ -1092,7 +1092,7 @@ static Index find_ordering /* return the number of garbage collections */ /* clear tag on pivot column */ Col [pivot_col].shared1.thickness = pivot_col_thickness ; - max_deg = (std::max)(max_deg, pivot_row_degree) ; + max_deg = numext::maxi(max_deg, pivot_row_degree) ; /* === Kill all rows used to construct pivot row ==================== */ @@ -1266,11 +1266,11 @@ static Index find_ordering /* return the number of garbage collections */ /* add set difference */ cur_score += row_mark - tag_mark ; /* integer overflow... */ - cur_score = (std::min)(cur_score, n_col) ; + cur_score = numext::mini(cur_score, n_col) ; } /* recompute the column's length */ - Col [col].length = (Index) (new_cp - &A [Col [col].start]) ; + Col [col].length = (IndexType) (new_cp - &A [Col [col].start]) ; /* === Further mass elimination ================================= */ @@ -1318,7 +1318,7 @@ static Index find_ordering /* return the number of garbage collections */ Col [col].shared4.hash_next = first_col ; /* save hash function in Col [col].shared3.hash */ - Col [col].shared3.hash = (Index) hash ; + Col [col].shared3.hash = (IndexType) hash ; COLAMD_ASSERT (COL_IS_ALIVE (col)) ; } } @@ -1379,7 +1379,7 @@ static Index find_ordering /* return the number of garbage collections */ cur_score -= Col [col].shared1.thickness ; /* make sure score is less or equal than the max score */ - cur_score = (std::min)(cur_score, max_score) ; + cur_score = numext::mini(cur_score, max_score) ; COLAMD_ASSERT (cur_score >= 0) ; /* store updated score */ @@ -1402,7 +1402,7 @@ static Index find_ordering /* return the number of garbage collections */ head [cur_score] = col ; /* see if this score is less than current min */ - min_score = (std::min)(min_score, cur_score) ; + min_score = numext::mini(min_score, cur_score) ; } @@ -1413,7 +1413,7 @@ static Index find_ordering /* return the number of garbage collections */ /* 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 = (Index) (new_rp - &A[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 */ @@ -1442,22 +1442,22 @@ static Index find_ordering /* return the number of garbage collections */ taken by this routine is O (n_col), that is, linear in the number of columns. Not user-callable. */ -template <typename Index> +template <typename IndexType> static inline void order_children ( /* === Parameters ======================================================= */ - Index n_col, /* number of columns of A */ - colamd_col<Index> Col [], /* of size n_col+1 */ - Index p [] /* p [0 ... n_col-1] is the column permutation*/ + IndexType n_col, /* number of columns of A */ + colamd_col<IndexType> Col [], /* of size n_col+1 */ + IndexType p [] /* p [0 ... n_col-1] is the column permutation*/ ) { /* === Local variables ================================================== */ - Index i ; /* loop counter for all columns */ - Index c ; /* column index */ - Index parent ; /* index of column's parent */ - Index order ; /* column's order */ + 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 ================================== */ @@ -1543,33 +1543,33 @@ static inline void order_children just been computed in the approximate degree computation. Not user-callable. */ -template <typename Index> +template <typename IndexType> static void detect_super_cols ( /* === Parameters ======================================================= */ - colamd_col<Index> Col [], /* of size n_col+1 */ - Index A [], /* row indices of A */ - Index head [], /* head of degree lists and hash buckets */ - Index row_start, /* pointer to set of columns to check */ - Index row_length /* number of columns to check */ + colamd_col<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 ================================================== */ - Index hash ; /* hash value for a column */ - Index *rp ; /* pointer to a row */ - Index c ; /* a column index */ - Index super_c ; /* column index of the column to absorb into */ - Index *cp1 ; /* column pointer for column super_c */ - Index *cp2 ; /* column pointer for column c */ - Index length ; /* length of column super_c */ - Index prev_c ; /* column preceding c in hash bucket */ - Index i ; /* loop counter */ - Index *rp_end ; /* pointer to the end of the row */ - Index col ; /* a column index in the row to check */ - Index head_column ; /* first column in hash bucket or degree list */ - Index first_col ; /* first column in hash bucket */ + 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 ================================== */ @@ -1694,27 +1694,27 @@ static void detect_super_cols itself linear in the number of nonzeros in the input matrix. Not user-callable. */ -template <typename Index> -static Index garbage_collection /* returns the new value of pfree */ +template <typename IndexType> +static IndexType garbage_collection /* returns the new value of pfree */ ( /* === Parameters ======================================================= */ - Index n_row, /* number of rows */ - Index n_col, /* number of columns */ - Colamd_Row<Index> Row [], /* row info */ - colamd_col<Index> Col [], /* column info */ - Index A [], /* A [0 ... Alen-1] holds the matrix */ - Index *pfree /* &A [0] ... pfree is in use */ + IndexType n_row, /* number of rows */ + IndexType n_col, /* number of columns */ + Colamd_Row<IndexType> Row [], /* row info */ + colamd_col<IndexType> Col [], /* column info */ + IndexType A [], /* A [0 ... Alen-1] holds the matrix */ + IndexType *pfree /* &A [0] ... pfree is in use */ ) { /* === Local variables ================================================== */ - Index *psrc ; /* source pointer */ - Index *pdest ; /* destination pointer */ - Index j ; /* counter */ - Index r ; /* a row index */ - Index c ; /* a column index */ - Index length ; /* length of a row or column */ + 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 =========================================== */ @@ -1727,7 +1727,7 @@ static Index garbage_collection /* returns the new value of pfree */ /* move and compact the column */ COLAMD_ASSERT (pdest <= psrc) ; - Col [c].start = (Index) (pdest - &A [0]) ; + Col [c].start = (IndexType) (pdest - &A [0]) ; length = Col [c].length ; for (j = 0 ; j < length ; j++) { @@ -1737,7 +1737,7 @@ static Index garbage_collection /* returns the new value of pfree */ *pdest++ = r ; } } - Col [c].length = (Index) (pdest - &A [Col [c].start]) ; + Col [c].length = (IndexType) (pdest - &A [Col [c].start]) ; } } @@ -1784,7 +1784,7 @@ static Index garbage_collection /* returns the new value of pfree */ /* move and compact the row */ COLAMD_ASSERT (pdest <= psrc) ; - Row [r].start = (Index) (pdest - &A [0]) ; + Row [r].start = (IndexType) (pdest - &A [0]) ; length = Row [r].length ; for (j = 0 ; j < length ; j++) { @@ -1794,7 +1794,7 @@ static Index garbage_collection /* returns the new value of pfree */ *pdest++ = c ; } } - Row [r].length = (Index) (pdest - &A [Row [r].start]) ; + Row [r].length = (IndexType) (pdest - &A [Row [r].start]) ; } } @@ -1803,7 +1803,7 @@ static Index garbage_collection /* returns the new value of pfree */ /* === Return the new value of pfree ==================================== */ - return ((Index) (pdest - &A [0])) ; + return ((IndexType) (pdest - &A [0])) ; } @@ -1815,18 +1815,18 @@ static Index garbage_collection /* returns the new value of pfree */ Clears the Row [].shared2.mark array, and returns the new tag_mark. Return value is the new tag_mark. Not user-callable. */ -template <typename Index> -static inline Index clear_mark /* return the new value for tag_mark */ +template <typename IndexType> +static inline IndexType clear_mark /* return the new value for tag_mark */ ( /* === Parameters ======================================================= */ - Index n_row, /* number of rows in A */ - Colamd_Row<Index> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ + IndexType n_row, /* number of rows in A */ + Colamd_Row<IndexType> Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ ) { /* === Local variables ================================================== */ - Index r ; + IndexType r ; for (r = 0 ; r < n_row ; r++) { diff --git a/eigen/Eigen/src/OrderingMethods/Ordering.h b/eigen/Eigen/src/OrderingMethods/Ordering.h index f3c31f9..7ea9b14 100644 --- a/eigen/Eigen/src/OrderingMethods/Ordering.h +++ b/eigen/Eigen/src/OrderingMethods/Ordering.h @@ -19,20 +19,21 @@ namespace internal { /** \internal * \ingroup OrderingMethods_Module - * \returns the symmetric pattern A^T+A from the input matrix A. + * \param[in] A the input non-symmetric matrix + * \param[out] symmat the symmetric pattern A^T+A from the input matrix \a A. * FIXME: The values should not be considered here */ template<typename MatrixType> -void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) +void ordering_helper_at_plus_a(const MatrixType& A, MatrixType& symmat) { MatrixType C; - C = mat.transpose(); // NOTE: Could be costly + C = A.transpose(); // NOTE: Could be costly for (int i = 0; i < C.rows(); i++) { for (typename MatrixType::InnerIterator it(C, i); it; ++it) it.valueRef() = 0.0; } - symmat = C + mat; + symmat = C + A; } } @@ -44,14 +45,14 @@ void ordering_helper_at_plus_a(const MatrixType& mat, MatrixType& symmat) * * Functor computing the \em approximate \em minimum \em degree ordering * If the matrix is not structurally symmetric, an ordering of A^T+A is computed - * \tparam Index The type of indices of the matrix + * \tparam StorageIndex The type of indices of the matrix * \sa COLAMDOrdering */ -template <typename Index> +template <typename StorageIndex> class AMDOrdering { public: - typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; + typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; /** Compute the permutation vector from a sparse matrix * This routine is much faster if the input matrix is column-major @@ -60,7 +61,7 @@ class AMDOrdering void operator()(const MatrixType& mat, PermutationType& perm) { // Compute the symmetric pattern - SparseMatrix<typename MatrixType::Scalar, ColMajor, Index> symm; + SparseMatrix<typename MatrixType::Scalar, ColMajor, StorageIndex> symm; internal::ordering_helper_at_plus_a(mat,symm); // Call the AMD routine @@ -72,7 +73,7 @@ class AMDOrdering template <typename SrcType, unsigned int SrcUpLo> void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm) { - SparseMatrix<typename SrcType::Scalar, ColMajor, Index> C; C = mat; + SparseMatrix<typename SrcType::Scalar, ColMajor, StorageIndex> C; C = mat; // Call the AMD routine // m_mat.prune(keep_diag()); //Remove the diagonal elements @@ -88,13 +89,13 @@ class AMDOrdering * Functor computing the natural ordering (identity) * * \note Returns an empty permutation matrix - * \tparam Index The type of indices of the matrix + * \tparam StorageIndex The type of indices of the matrix */ -template <typename Index> +template <typename StorageIndex> class NaturalOrdering { public: - typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; + typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; /** Compute the permutation vector from a column-major sparse matrix */ template <typename MatrixType> @@ -108,15 +109,17 @@ class NaturalOrdering /** \ingroup OrderingMethods_Module * \class COLAMDOrdering * + * \tparam StorageIndex The type of indices of the matrix + * * Functor computing the \em column \em approximate \em minimum \em degree ordering * The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()). */ -template<typename Index> +template<typename StorageIndex> class COLAMDOrdering { public: - typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; - typedef Matrix<Index, Dynamic, 1> IndexVector; + typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; + typedef Matrix<StorageIndex, Dynamic, 1> IndexVector; /** Compute the permutation vector \a perm form the sparse matrix \a mat * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). @@ -126,26 +129,26 @@ class COLAMDOrdering { eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering"); - Index m = mat.rows(); - Index n = mat.cols(); - Index nnz = mat.nonZeros(); + StorageIndex m = StorageIndex(mat.rows()); + StorageIndex n = StorageIndex(mat.cols()); + StorageIndex nnz = StorageIndex(mat.nonZeros()); // Get the recommended value of Alen to be used by colamd - Index Alen = internal::colamd_recommended(nnz, m, n); + StorageIndex Alen = internal::colamd_recommended(nnz, m, n); // Set the default parameters double knobs [COLAMD_KNOBS]; - Index stats [COLAMD_STATS]; + StorageIndex stats [COLAMD_STATS]; internal::colamd_set_defaults(knobs); IndexVector p(n+1), A(Alen); - for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; - for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; + for(StorageIndex i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; + for(StorageIndex i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; // Call Colamd routine to compute the ordering - Index info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + StorageIndex info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); EIGEN_UNUSED_VARIABLE(info); eigen_assert( info && "COLAMD failed " ); perm.resize(n); - for (Index i = 0; i < n; i++) perm.indices()(p(i)) = i; + for (StorageIndex i = 0; i < n; i++) perm.indices()(p(i)) = i; } }; |