add SparseLU in sparse bench
diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h index 0af137d..686c0f9 100644 --- a/Eigen/src/OrderingMethods/Eigen_Colamd.h +++ b/Eigen/src/OrderingMethods/Eigen_Colamd.h
@@ -7,7 +7,7 @@ // 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 eigen_colamd/symamd library. The copyright is below +// 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 @@ -39,18 +39,19 @@ // // Availability: // -// The eigen_colamd/symamd library is available at +// The colamd/symamd library is available at // -// http://www.cise.ufl.edu/research/sparse/eigen_colamd/ +// http://www.cise.ufl.edu/research/sparse/colamd/ -// This is the http://www.cise.ufl.edu/research/sparse/eigen_colamd/eigen_colamd.h -// file. It is required by the eigen_colamd.c, colamdmex.c, and symamdmex.c +// This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h +// file. It is required by the colamd.c, colamdmex.c, and symamdmex.c // files, and by any C code that calls the routines whose prototypes are -// listed below, or that uses the eigen_colamd/symamd definitions listed below. +// listed below, or that uses the colamd/symamd definitions listed below. #ifndef EIGEN_COLAMD_H #define EIGEN_COLAMD_H +namespace internal { /* Ensure that debugging is turned off: */ #ifndef COLAMD_NDEBUG #define COLAMD_NDEBUG @@ -61,42 +62,42 @@ /* ========================================================================== */ /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ -#define EIGEN_COLAMD_KNOBS 20 +#define COLAMD_KNOBS 20 /* number of output statistics. Only stats [0..6] are currently used. */ -#define EIGEN_COLAMD_STATS 20 +#define COLAMD_STATS 20 /* knobs [0] and stats [0]: dense row knob and output statistic. */ -#define EIGEN_COLAMD_DENSE_ROW 0 +#define COLAMD_DENSE_ROW 0 /* knobs [1] and stats [1]: dense column knob and output statistic. */ -#define EIGEN_COLAMD_DENSE_COL 1 +#define COLAMD_DENSE_COL 1 /* stats [2]: memory defragmentation count output statistic */ -#define EIGEN_COLAMD_DEFRAG_COUNT 2 +#define COLAMD_DEFRAG_COUNT 2 -/* stats [3]: eigen_colamd status: zero OK, > 0 warning or notice, < 0 error */ -#define EIGEN_COLAMD_STATUS 3 +/* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ +#define COLAMD_STATUS 3 /* stats [4..6]: error info, or info on jumbled columns */ -#define EIGEN_COLAMD_INFO1 4 -#define EIGEN_COLAMD_INFO2 5 -#define EIGEN_COLAMD_INFO3 6 +#define COLAMD_INFO1 4 +#define COLAMD_INFO2 5 +#define COLAMD_INFO3 6 /* error codes returned in stats [3]: */ -#define EIGEN_COLAMD_OK (0) -#define EIGEN_COLAMD_OK_BUT_JUMBLED (1) -#define EIGEN_COLAMD_ERROR_A_not_present (-1) -#define EIGEN_COLAMD_ERROR_p_not_present (-2) -#define EIGEN_COLAMD_ERROR_nrow_negative (-3) -#define EIGEN_COLAMD_ERROR_ncol_negative (-4) -#define EIGEN_COLAMD_ERROR_nnz_negative (-5) -#define EIGEN_COLAMD_ERROR_p0_nonzero (-6) -#define EIGEN_COLAMD_ERROR_A_too_small (-7) -#define EIGEN_COLAMD_ERROR_col_length_negative (-8) -#define EIGEN_COLAMD_ERROR_row_index_out_of_bounds (-9) -#define EIGEN_COLAMD_ERROR_out_of_memory (-10) -#define EIGEN_COLAMD_ERROR_internal_error (-999) +#define COLAMD_OK (0) +#define COLAMD_OK_BUT_JUMBLED (1) +#define COLAMD_ERROR_A_not_present (-1) +#define COLAMD_ERROR_p_not_present (-2) +#define COLAMD_ERROR_nrow_negative (-3) +#define COLAMD_ERROR_ncol_negative (-4) +#define COLAMD_ERROR_nnz_negative (-5) +#define COLAMD_ERROR_p0_nonzero (-6) +#define COLAMD_ERROR_A_too_small (-7) +#define COLAMD_ERROR_col_length_negative (-8) +#define COLAMD_ERROR_row_index_out_of_bounds (-9) +#define COLAMD_ERROR_out_of_memory (-10) +#define COLAMD_ERROR_internal_error (-999) /* ========================================================================== */ /* === Definitions ========================================================== */ @@ -105,30 +106,30 @@ #define COLAMD_MAX(a,b) (((a) > (b)) ? (a) : (b)) #define COLAMD_MIN(a,b) (((a) < (b)) ? (a) : (b)) -#define EIGEN_ONES_COMPLEMENT(r) (-(r)-1) +#define ONES_COMPLEMENT(r) (-(r)-1) /* -------------------------------------------------------------------------- */ -#define EIGEN_COLAMD_EMPTY (-1) +#define COLAMD_EMPTY (-1) /* Row and column status */ -#define EIGEN_ALIVE (0) -#define EIGEN_DEAD (-1) +#define ALIVE (0) +#define DEAD (-1) /* Column status */ -#define EIGEN_DEAD_PRINCIPAL (-1) -#define EIGEN_DEAD_NON_PRINCIPAL (-2) +#define DEAD_PRINCIPAL (-1) +#define DEAD_NON_PRINCIPAL (-2) /* Macros for row and column status update and checking. */ -#define EIGEN_ROW_IS_DEAD(r) EIGEN_ROW_IS_MARKED_DEAD (Row[r].shared2.mark) -#define EIGEN_ROW_IS_MARKED_DEAD(row_mark) (row_mark < EIGEN_ALIVE) -#define EIGEN_ROW_IS_ALIVE(r) (Row [r].shared2.mark >= EIGEN_ALIVE) -#define EIGEN_COL_IS_DEAD(c) (Col [c].start < EIGEN_ALIVE) -#define EIGEN_COL_IS_ALIVE(c) (Col [c].start >= EIGEN_ALIVE) -#define EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == EIGEN_DEAD_PRINCIPAL) -#define EIGEN_KILL_ROW(r) { Row [r].shared2.mark = EIGEN_DEAD ; } -#define EIGEN_KILL_PRINCIPAL_COL(c) { Col [c].start = EIGEN_DEAD_PRINCIPAL ; } -#define EIGEN_KILL_NON_PRINCIPAL_COL(c) { Col [c].start = EIGEN_DEAD_NON_PRINCIPAL ; } +#define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) +#define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) +#define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) +#define COL_IS_DEAD(c) (Col [c].start < ALIVE) +#define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) +#define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) +#define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } +#define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } +#define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } /* ========================================================================== */ /* === Colamd reporting mechanism =========================================== */ @@ -146,7 +147,7 @@ /* Use printf in standard C environment, for debugging and statistics output. */ /* Output is generated only if debugging is enabled at compile time, or if */ -/* the caller explicitly calls eigen_colamd_report or symamd_report. */ +/* the caller explicitly calls colamd_report or symamd_report. */ #define PRINTF printf /* In C, matrices are 0-based and indices are reported as such in *_report */ @@ -155,9 +156,9 @@ #endif /* MATLAB_MEX_FILE */ // == Row and Column structures == -typedef struct EIGEN_Colamd_Col_struct +typedef struct colamd_col_struct { - int start ; /* index for A of first row in this column, or EIGEN_DEAD */ + int start ; /* index for A of first row in this column, or DEAD */ /* if column is dead */ int length ; /* number of rows in this column */ union @@ -186,16 +187,16 @@ int hash_next ; /* next column, if col is in a hash list */ } shared4 ; -} EIGEN_Colamd_Col ; +} colamd_col ; -typedef struct EIGEN_Colamd_Row_struct +typedef struct Colamd_Row_struct { int start ; /* index for A of first col in this row */ int length ; /* number of principal columns in this row */ union { int degree ; /* number of principal & non-principal columns in row */ - int p ; /* used as a row pointer in eigen_init_rows_cols () */ + int p ; /* used as a row pointer in init_rows_cols () */ } shared1 ; union { @@ -203,19 +204,19 @@ int first_column ;/* first column in row (used in garbage collection) */ } shared2 ; -} EIGEN_Colamd_Row ; +} Colamd_Row ; /* ========================================================================== */ /* === Colamd recommended memory size ======================================= */ /* ========================================================================== */ /* - The recommended length Alen of the array A passed to eigen_colamd is given by - the EIGEN_COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any + 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. EIGEN_COLAMD_C (n_col) + EIGEN_COLAMD_R (n_row) space is + 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 - eigen_colamd. An additional n_col space is the minimal amount of "elbow room", + 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. @@ -224,120 +225,41 @@ gcc -pedantic warning messages. */ -#define EIGEN_COLAMD_C(n_col) ((int) (((n_col) + 1) * sizeof (EIGEN_Colamd_Col) / sizeof (int))) -#define EIGEN_COLAMD_R(n_row) ((int) (((n_row) + 1) * sizeof (EIGEN_Colamd_Row) / sizeof (int))) +inline int colamd_c(int n_col) +{ return int( ((n_col) + 1) * sizeof (colamd_col) / sizeof (int) ) ; } -#define EIGEN_COLAMD_RECOMMENDED(nnz, n_row, n_col) \ -( \ -((nnz) < 0 || (n_row) < 0 || (n_col) < 0) \ -? \ - (-1) \ -: \ - (2 * (nnz) + EIGEN_COLAMD_C (n_col) + EIGEN_COLAMD_R (n_row) + (n_col) + ((nnz) / 5)) \ -) +inline int colamd_r(int n_row) +{ return int(((n_row) + 1) * sizeof (Colamd_Row) / sizeof (int)); } // Various routines -int eigen_colamd_recommended (int nnz, int n_row, int n_col) ; +inline int colamd_recommended (int nnz, int n_row, int n_col) ; -void eigen_colamd_set_defaults (double knobs [EIGEN_COLAMD_KNOBS]) ; +static inline void colamd_set_defaults (double knobs [COLAMD_KNOBS]) ; -bool eigen_colamd (int n_row, int n_col, int Alen, int A [], int p [], double knobs[EIGEN_COLAMD_KNOBS], int stats [EIGEN_COLAMD_STATS]) ; +static bool colamd (int n_row, int n_col, int Alen, int A [], int p [], double knobs[COLAMD_KNOBS], int stats [COLAMD_STATS]) ; -void eigen_colamd_report (int stats [EIGEN_COLAMD_STATS]); +static inline void colamd_report (int stats [COLAMD_STATS]); -int eigen_init_rows_cols (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col col [], int A [], int p [], int stats[EIGEN_COLAMD_STATS] ); +static int init_rows_cols (int n_row, int n_col, Colamd_Row Row [], colamd_col col [], int A [], int p [], int stats[COLAMD_STATS] ); -void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); +static void init_scoring (int n_row, int n_col, Colamd_Row Row [], colamd_col Col [], int A [], int head [], double knobs[COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); -int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], int n_col2, int max_deg, int pfree); +static int find_ordering (int n_row, int n_col, int Alen, Colamd_Row Row [], colamd_col Col [], int A [], int head [], int n_col2, int max_deg, int pfree); -void eigen_order_children (int n_col, EIGEN_Colamd_Col Col [], int p []); +static void order_children (int n_col, colamd_col Col [], int p []); -void eigen_detect_super_cols ( -#ifndef COLAMD_NDEBUG - int n_col, - EIGEN_Colamd_Row Row [], -#endif /* COLAMD_NDEBUG */ - EIGEN_Colamd_Col Col [], +static void detect_super_cols ( + colamd_col Col [], int A [], int head [], int row_start, int row_length ) ; - int eigen_garbage_collection (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int *pfree) ; +static int garbage_collection (int n_row, int n_col, Colamd_Row Row [], colamd_col Col [], int A [], int *pfree) ; - int eigen_clear_mark (int n_row, EIGEN_Colamd_Row Row [] ) ; +static inline int clear_mark (int n_row, Colamd_Row Row [] ) ; - void eigen_print_report (char *method, int stats [EIGEN_COLAMD_STATS]) ; - -/* ========================================================================== */ -/* === Debugging prototypes and definitions ================================= */ -/* ========================================================================== */ - -#ifndef COLAMD_NDEBUG - -/* colamd_debug is the *ONLY* global variable, and is only */ -/* present when debugging */ - - int colamd_debug ; /* debug print level */ - -#define COLAMD_DEBUG0(params) { (void) PRINTF params ; } -#define COLAMD_DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; } -#define COLAMD_DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; } -#define COLAMD_DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; } -#define COLAMD_DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; } - -#ifdef MATLAB_MEX_FILE -#define COLAMD_ASSERT(expression) (mxAssert ((expression), "")) -#else -#define COLAMD_ASSERT(expression) (assert (expression)) -#endif /* MATLAB_MEX_FILE */ - - void eigen_colamd_get_debug /* gets the debug print level from getenv */ -( - char *method -) ; - - void eigen_debug_deg_lists -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int head [], - int min_score, - int should, - int max_deg -) ; - - void eigen_debug_mark -( - int n_row, - EIGEN_Colamd_Row Row [], - int tag_mark, - int max_mark -) ; - - void eigen_debug_matrix -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [] -) ; - - void eigen_debug_structures -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [], - int n_col2 -) ; - -#else /* COLAMD_NDEBUG */ +static void print_report (const char *method, int stats [COLAMD_STATS]) ; /* === No debugging ========================================================= */ @@ -349,51 +271,50 @@ #define COLAMD_ASSERT(expression) ((void) 0) -#endif /* COLAMD_NDEBUG */ - - /** * \brief Returns the recommended value of Alen * - * Returns recommended value of Alen for use by eigen_colamd. + * 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 EIGEN_COLAMD_RECOMMENDED. + * 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 eigen_colamd + * \return recommended value of Alen for use by colamd */ -int eigen_colamd_recommended ( int nnz, int n_row, int n_col) +inline int colamd_recommended ( int nnz, int n_row, int n_col) { - - return (EIGEN_COLAMD_RECOMMENDED (nnz, n_row, 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 [EIGEN_COLAMD_DENSE_ROW] * n_col) + * Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) * entries are removed prior to ordering. Columns with more than - * (knobs [EIGEN_COLAMD_DENSE_COL] * n_row) entries are removed prior to + * (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to * ordering, and placed last in the output column ordering. * - * EIGEN_COLAMD_DENSE_ROW and EIGEN_COLAMD_DENSE_COL are defined as 0 and 1, - * respectively, in eigen_colamd.h. Default values of these two knobs + * COLAMD_DENSE_ROW and COLAMD_DENSE_COL 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 - * eigen_colamd_set_defaults, so that the code that calls eigen_colamd will + * colamd_set_defaults, so that the code that calls colamd will * not need to change, assuming that you either use - * eigen_colamd_set_defaults, or pass a (double *) NULL pointer as the - * knobs array to eigen_colamd or symamd. + * colamd_set_defaults, or pass a (double *) NULL pointer as the + * knobs array to colamd or symamd. * - * \param knobs parameter settings for eigen_colamd + * \param knobs parameter settings for colamd */ -void eigen_colamd_set_defaults(double knobs[EIGEN_COLAMD_KNOBS]) +static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS]) { /* === Local variables ================================================== */ @@ -403,12 +324,12 @@ { return ; /* no knobs to initialize */ } - for (i = 0 ; i < EIGEN_COLAMD_KNOBS ; i++) + for (i = 0 ; i < COLAMD_KNOBS ; i++) { knobs [i] = 0 ; } - knobs [EIGEN_COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ - knobs [EIGEN_COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ + knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ + knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ } /** @@ -425,10 +346,10 @@ * \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 eigen_colamd - * \param stats eigen_colamd output statistics and error codes + * \param knobs parameter settings for colamd + * \param stats colamd output statistics and error codes */ -bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[EIGEN_COLAMD_KNOBS], int stats[EIGEN_COLAMD_STATS]) +static bool colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[COLAMD_KNOBS], int stats[COLAMD_STATS]) { /* === Local variables ================================================== */ @@ -437,77 +358,74 @@ int Row_size ; /* size of Row [], in integers */ int Col_size ; /* size of Col [], in integers */ int need ; /* minimum required length of A */ - EIGEN_Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ - EIGEN_Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ + Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ + colamd_col *Col ; /* pointer into A of Col [0..n_col] array */ int n_col2 ; /* number of non-dense, non-empty columns */ int n_row2 ; /* number of non-dense, non-empty rows */ int ngarbage ; /* number of garbage collections performed */ int max_deg ; /* maximum row degree */ - double default_knobs [EIGEN_COLAMD_KNOBS] ; /* default knobs array */ + double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ -#ifndef COLAMD_NDEBUG - eigen_colamd_get_debug ("eigen_colamd") ; -#endif /* COLAMD_NDEBUG */ /* === Check the input arguments ======================================== */ if (!stats) { - COLAMD_DEBUG0 (("eigen_colamd: stats not present\n")) ; + COLAMD_DEBUG0 (("colamd: stats not present\n")) ; return (false) ; } - for (i = 0 ; i < EIGEN_COLAMD_STATS ; i++) + for (i = 0 ; i < COLAMD_STATS ; i++) { stats [i] = 0 ; } - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_OK ; - stats [EIGEN_COLAMD_INFO1] = -1 ; - stats [EIGEN_COLAMD_INFO2] = -1 ; + stats [COLAMD_STATUS] = COLAMD_OK ; + stats [COLAMD_INFO1] = -1 ; + stats [COLAMD_INFO2] = -1 ; if (!A) /* A is not present */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_A_not_present ; - COLAMD_DEBUG0 (("eigen_colamd: A not present\n")) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; + COLAMD_DEBUG0 (("colamd: A not present\n")) ; return (false) ; } if (!p) /* p is not present */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_p_not_present ; - COLAMD_DEBUG0 (("eigen_colamd: p not present\n")) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; + COLAMD_DEBUG0 (("colamd: p not present\n")) ; return (false) ; } if (n_row < 0) /* n_row must be >= 0 */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_nrow_negative ; - stats [EIGEN_COLAMD_INFO1] = n_row ; - COLAMD_DEBUG0 (("eigen_colamd: nrow negative %d\n", n_row)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; + 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 [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_ncol_negative ; - stats [EIGEN_COLAMD_INFO1] = n_col ; - COLAMD_DEBUG0 (("eigen_colamd: ncol negative %d\n", n_col)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; + 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 [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_nnz_negative ; - stats [EIGEN_COLAMD_INFO1] = nnz ; - COLAMD_DEBUG0 (("eigen_colamd: number of entries negative %d\n", nnz)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; + stats [COLAMD_INFO1] = nnz ; + COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; return (false) ; } if (p [0] != 0) { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_p0_nonzero ; - stats [EIGEN_COLAMD_INFO1] = p [0] ; - COLAMD_DEBUG0 (("eigen_colamd: p[0] not zero %d\n", p [0])) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; + stats [COLAMD_INFO1] = p [0] ; + COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; return (false) ; } @@ -515,72 +433,73 @@ if (!knobs) { - eigen_colamd_set_defaults (default_knobs) ; + colamd_set_defaults (default_knobs) ; knobs = default_knobs ; } /* === Allocate the Row and Col arrays from array A ===================== */ - Col_size = EIGEN_COLAMD_C (n_col) ; - Row_size = EIGEN_COLAMD_R (n_row) ; + 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 [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_A_too_small ; - stats [EIGEN_COLAMD_INFO1] = need ; - stats [EIGEN_COLAMD_INFO2] = Alen ; - COLAMD_DEBUG0 (("eigen_colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); + stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; + 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 = (EIGEN_Colamd_Col *) &A [Alen] ; - Row = (EIGEN_Colamd_Row *) &A [Alen + Col_size] ; + Col = (colamd_col *) &A [Alen] ; + Row = (Colamd_Row *) &A [Alen + Col_size] ; /* === Construct the row and column data structures ===================== */ - if (!eigen_init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) + if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) { /* input matrix is invalid */ - COLAMD_DEBUG0 (("eigen_colamd: Matrix invalid\n")) ; + COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ; return (false) ; } /* === Initialize scores, kill dense rows/columns ======================= */ - eigen_init_scoring (n_row, n_col, Row, Col, A, p, knobs, + init_scoring (n_row, n_col, Row, Col, A, p, knobs, &n_row2, &n_col2, &max_deg) ; /* === Order the supercolumns =========================================== */ - ngarbage = eigen_find_ordering (n_row, n_col, Alen, Row, Col, A, p, + ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, n_col2, max_deg, 2*nnz) ; /* === Order the non-principal columns ================================== */ - eigen_order_children (n_col, Col, p) ; + order_children (n_col, Col, p) ; /* === Return statistics in stats ======================================= */ - stats [EIGEN_COLAMD_DENSE_ROW] = n_row - n_row2 ; - stats [EIGEN_COLAMD_DENSE_COL] = n_col - n_col2 ; - stats [EIGEN_COLAMD_DEFRAG_COUNT] = ngarbage ; - COLAMD_DEBUG0 (("eigen_colamd: done.\n")) ; + stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; + stats [COLAMD_DENSE_COL] = n_col - n_col2 ; + stats [COLAMD_DEFRAG_COUNT] = ngarbage ; + COLAMD_DEBUG0 (("colamd: done.\n")) ; return (true) ; } /* ========================================================================== */ -/* === eigen_colamd_report ======================================================== */ +/* === colamd_report ======================================================== */ /* ========================================================================== */ - void eigen_colamd_report + static inline void colamd_report ( - int stats [EIGEN_COLAMD_STATS] + int stats [COLAMD_STATS] ) { - eigen_print_report ("eigen_colamd", stats) ; + const char *method = "colamd"; + print_report (method, stats) ; } @@ -592,7 +511,7 @@ /* ========================================================================== */ -/* === eigen_init_rows_cols ======================================================= */ +/* === init_rows_cols ======================================================= */ /* ========================================================================== */ /* @@ -604,17 +523,17 @@ true otherwise. Not user-callable. */ - int eigen_init_rows_cols /* returns true if OK, or false otherwise */ + static int init_rows_cols /* returns true if OK, or false otherwise */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* row indices of A, of size Alen */ int p [], /* pointers to columns in A, of size n_col+1 */ - int stats [EIGEN_COLAMD_STATS] /* eigen_colamd statistics */ + int stats [COLAMD_STATS] /* colamd statistics */ ) { /* === Local variables ================================================== */ @@ -637,24 +556,24 @@ if (Col [col].length < 0) { /* column pointers must be non-decreasing */ - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_col_length_negative ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = Col [col].length ; - COLAMD_DEBUG0 (("eigen_colamd: col %d length %d < 0\n", col, Col [col].length)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; + 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 = EIGEN_COLAMD_EMPTY ; - Col [col].shared4.degree_next = EIGEN_COLAMD_EMPTY ; + Col [col].shared3.prev = COLAMD_EMPTY ; + Col [col].shared4.degree_next = COLAMD_EMPTY ; } /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ /* === Scan columns, compute row degrees, and check row indices ========= */ - stats [EIGEN_COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ + stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ for (row = 0 ; row < n_row ; row++) { @@ -676,11 +595,11 @@ /* make sure row indices within range */ if (row < 0 || row >= n_row) { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_row_index_out_of_bounds ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = row ; - stats [EIGEN_COLAMD_INFO3] = n_row ; - COLAMD_DEBUG0 (("eigen_colamd: row %d col %d out of bounds\n", row, col)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; + 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) ; } @@ -688,11 +607,11 @@ { /* row index are unsorted or repeated (or both), thus col */ /* is jumbled. This is a notice, not an error condition. */ - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_OK_BUT_JUMBLED ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = row ; - (stats [EIGEN_COLAMD_INFO3]) ++ ; - COLAMD_DEBUG1 (("eigen_colamd: row %d col %d unsorted/duplicate\n",row,col)); + stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; + 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) @@ -729,7 +648,7 @@ /* === Create row form ================================================== */ - if (stats [EIGEN_COLAMD_STATUS] == EIGEN_COLAMD_OK_BUT_JUMBLED) + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { /* if cols jumbled, watch for repeated row indices */ for (col = 0 ; col < n_col ; col++) @@ -771,31 +690,10 @@ /* === See if we need to re-create columns ============================== */ - if (stats [EIGEN_COLAMD_STATUS] == EIGEN_COLAMD_OK_BUT_JUMBLED) + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { - COLAMD_DEBUG0 (("eigen_colamd: reconstructing column form, matrix jumbled\n")) ; + COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; -#ifndef COLAMD_NDEBUG - /* make sure column lengths are correct */ - for (col = 0 ; col < n_col ; col++) - { - p [col] = Col [col].length ; - } - for (row = 0 ; row < n_row ; row++) - { - rp = &A [Row [row].start] ; - rp_end = rp + Row [row].length ; - while (rp < rp_end) - { - p [*rp++]-- ; - } - } - for (col = 0 ; col < n_col ; col++) - { - COLAMD_ASSERT (p [col] == 0) ; - } - /* now p is all zero (different than when debugging is turned off) */ -#endif /* COLAMD_NDEBUG */ /* === Compute col pointers ========================================= */ @@ -833,7 +731,7 @@ /* ========================================================================== */ -/* === eigen_init_scoring ========================================================= */ +/* === init_scoring ========================================================= */ /* ========================================================================== */ /* @@ -841,17 +739,17 @@ each column, and places all columns in the degree lists. Not user-callable. */ - void eigen_init_scoring +static void init_scoring ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ - double knobs [EIGEN_COLAMD_KNOBS],/* parameters */ + double knobs [COLAMD_KNOBS],/* parameters */ int *p_n_row2, /* number of non-dense, non-empty rows */ int *p_n_col2, /* number of non-dense, non-empty columns */ int *p_max_deg /* maximum row degree */ @@ -875,15 +773,12 @@ int max_deg ; /* maximum row degree */ int next_col ; /* Used to add to degree list.*/ -#ifndef COLAMD_NDEBUG - int debug_count ; /* debug only. */ -#endif /* COLAMD_NDEBUG */ /* === Extract knobs ==================================================== */ - dense_row_count = COLAMD_MAX (0, COLAMD_MIN (knobs [EIGEN_COLAMD_DENSE_ROW] * n_col, n_col)) ; - dense_col_count = COLAMD_MAX (0, COLAMD_MIN (knobs [EIGEN_COLAMD_DENSE_COL] * n_row, n_row)) ; - COLAMD_DEBUG1 (("eigen_colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; + dense_row_count = COLAMD_MAX (0, COLAMD_MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ; + dense_col_count = COLAMD_MAX (0, COLAMD_MIN (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 ; n_row2 = n_row ; @@ -899,10 +794,10 @@ { /* this is a empty column, kill and order it last */ Col [c].shared2.order = --n_col2 ; - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } } - COLAMD_DEBUG1 (("eigen_colamd: null columns killed: %d\n", n_col - n_col2)) ; + COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense columns =============================================== */ @@ -910,7 +805,7 @@ for (c = n_col-1 ; c >= 0 ; c--) { /* skip any dead columns */ - if (EIGEN_COL_IS_DEAD (c)) + if (COL_IS_DEAD (c)) { continue ; } @@ -926,10 +821,10 @@ { Row [*cp++].shared1.degree-- ; } - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; + COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense and empty rows ======================================== */ @@ -940,7 +835,7 @@ if (deg > dense_row_count || deg == 0) { /* kill a dense or empty row */ - EIGEN_KILL_ROW (r) ; + KILL_ROW (r) ; --n_row2 ; } else @@ -949,7 +844,7 @@ max_deg = COLAMD_MAX (max_deg, deg) ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; + COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; /* === Compute initial column scores ==================================== */ @@ -962,7 +857,7 @@ for (c = n_col-1 ; c >= 0 ; c--) { /* skip dead column */ - if (EIGEN_COL_IS_DEAD (c)) + if (COL_IS_DEAD (c)) { continue ; } @@ -975,7 +870,7 @@ /* get a row */ row = *cp++ ; /* skip if dead */ - if (EIGEN_ROW_IS_DEAD (row)) + if (ROW_IS_DEAD (row)) { continue ; } @@ -994,7 +889,7 @@ /* and have already been killed) */ COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ; Col [c].shared2.order = --n_col2 ; - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } else { @@ -1005,7 +900,7 @@ Col [c].shared2.score = score ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense, null, and newly-null columns killed: %d\n", + 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 */ @@ -1013,20 +908,13 @@ /* yet). Rows may contain dead columns, but all live rows contain at */ /* least one live column. */ -#ifndef COLAMD_NDEBUG - eigen_debug_structures (n_row, n_col, Row, Col, A, n_col2) ; -#endif /* COLAMD_NDEBUG */ - /* === Initialize degree lists ========================================== */ -#ifndef COLAMD_NDEBUG - debug_count = 0 ; -#endif /* COLAMD_NDEBUG */ /* clear the hash buckets */ for (c = 0 ; c <= n_col ; c++) { - head [c] = EIGEN_COLAMD_EMPTY ; + head [c] = COLAMD_EMPTY ; } min_score = n_col ; /* place in reverse order, so low column indices are at the front */ @@ -1034,7 +922,7 @@ for (c = n_col-1 ; c >= 0 ; c--) { /* only add principal columns to degree lists */ - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n", c, Col [c].shared2.score, min_score, n_col)) ; @@ -1047,16 +935,16 @@ COLAMD_ASSERT (min_score <= n_col) ; COLAMD_ASSERT (score >= 0) ; COLAMD_ASSERT (score <= n_col) ; - COLAMD_ASSERT (head [score] >= EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ; /* now add this column to dList at proper score location */ next_col = head [score] ; - Col [c].shared3.prev = EIGEN_COLAMD_EMPTY ; + Col [c].shared3.prev = COLAMD_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 != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = c ; } @@ -1065,19 +953,10 @@ /* see if this score is less than current min */ min_score = COLAMD_MIN (min_score, score) ; -#ifndef COLAMD_NDEBUG - debug_count++ ; -#endif /* COLAMD_NDEBUG */ } } -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG1 (("eigen_colamd: Live cols %d out of %d, non-princ: %d\n", - debug_count, n_col, n_col-debug_count)) ; - COLAMD_ASSERT (debug_count == n_col2) ; - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; -#endif /* COLAMD_NDEBUG */ /* === Return number of remaining columns, and max row degree =========== */ @@ -1088,7 +967,7 @@ /* ========================================================================== */ -/* === eigen_find_ordering ======================================================== */ +/* === find_ordering ======================================================== */ /* ========================================================================== */ /* @@ -1097,15 +976,15 @@ degree ordering method. Not user-callable. */ - int eigen_find_ordering /* return the number of garbage collections */ +static int find_ordering /* return the number of garbage collections */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ int Alen, /* size of A, 2*nnz + n_col or larger */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ int n_col2, /* Remaining columns to order */ @@ -1147,55 +1026,29 @@ int next_col ; /* Used by Dlist operations. */ int ngarbage ; /* number of garbage collections performed */ -#ifndef COLAMD_NDEBUG - int debug_d ; /* debug loop counter */ - int debug_step = 0 ; /* debug loop counter */ -#endif /* COLAMD_NDEBUG */ /* === Initialization and clear mark ==================================== */ max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */ - tag_mark = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; min_score = 0 ; ngarbage = 0 ; - COLAMD_DEBUG1 (("eigen_colamd: Ordering, n_col2=%d\n", n_col2)) ; + COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; /* === Order the columns ================================================ */ for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) { -#ifndef COLAMD_NDEBUG - if (debug_step % 100 == 0) - { - COLAMD_DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; - } - else - { - COLAMD_DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; - } - debug_step++ ; - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k, max_deg) ; - eigen_debug_matrix (n_row, n_col, Row, Col, A) ; -#endif /* COLAMD_NDEBUG */ - /* === 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] >= EIGEN_COLAMD_EMPTY) ; - -#ifndef COLAMD_NDEBUG - for (debug_d = 0 ; debug_d < min_score ; debug_d++) - { - COLAMD_ASSERT (head [debug_d] == EIGEN_COLAMD_EMPTY) ; - } -#endif /* COLAMD_NDEBUG */ + COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ; /* get pivot column from head of minimum degree list */ - while (head [min_score] == EIGEN_COLAMD_EMPTY && min_score < n_col) + while (head [min_score] == COLAMD_EMPTY && min_score < n_col) { min_score++ ; } @@ -1203,12 +1056,12 @@ 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 != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { - Col [next_col].shared3.prev = EIGEN_COLAMD_EMPTY ; + Col [next_col].shared3.prev = COLAMD_EMPTY ; } - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (pivot_col)) ; + COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ; COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ; /* remember score for defrag check */ @@ -1227,16 +1080,13 @@ needed_memory = COLAMD_MIN (pivot_col_score, n_col - k) ; if (pfree + needed_memory >= Alen) { - pfree = eigen_garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; + pfree = 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 = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; -#ifndef COLAMD_NDEBUG - eigen_debug_matrix (n_row, n_col, Row, Col, A) ; -#endif /* COLAMD_NDEBUG */ } /* === Compute pivot row pattern ==================================== */ @@ -1258,9 +1108,9 @@ { /* get a row */ row = *cp++ ; - COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", EIGEN_ROW_IS_ALIVE (row), row)) ; + COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; /* skip if row is dead */ - if (EIGEN_ROW_IS_DEAD (row)) + if (ROW_IS_DEAD (row)) { continue ; } @@ -1272,7 +1122,7 @@ col = *rp++ ; /* add the column, if alive and untagged */ col_thickness = Col [col].shared1.thickness ; - if (col_thickness > 0 && EIGEN_COL_IS_ALIVE (col)) + if (col_thickness > 0 && COL_IS_ALIVE (col)) { /* tag column in pivot row */ Col [col].shared1.thickness = -col_thickness ; @@ -1288,10 +1138,6 @@ Col [pivot_col].shared1.thickness = pivot_col_thickness ; max_deg = COLAMD_MAX (max_deg, pivot_row_degree) ; -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG3 (("check2\n")) ; - eigen_debug_mark (n_row, Row, tag_mark, max_mark) ; -#endif /* COLAMD_NDEBUG */ /* === Kill all rows used to construct pivot row ==================== */ @@ -1303,7 +1149,7 @@ /* may be killing an already dead row */ row = *cp++ ; COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ; - EIGEN_KILL_ROW (row) ; + KILL_ROW (row) ; } /* === Select a row index to use as the new pivot row =============== */ @@ -1318,7 +1164,7 @@ else { /* there is no pivot row, since it is of zero length */ - pivot_row = EIGEN_COLAMD_EMPTY ; + pivot_row = COLAMD_EMPTY ; COLAMD_ASSERT (pivot_row_length == 0) ; } COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; @@ -1340,7 +1186,7 @@ /* 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 eigen_colamd. */ + /* column on input to colamd. */ /* === Compute set differences ====================================== */ @@ -1355,7 +1201,7 @@ while (rp < rp_end) { col = *rp++ ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col) && col != pivot_col) ; + COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; COLAMD_DEBUG3 (("Col: %d\n", col)) ; /* clear tags used to construct pivot row pattern */ @@ -1370,8 +1216,8 @@ next_col = Col [col].shared4.degree_next ; COLAMD_ASSERT (cur_score >= 0) ; COLAMD_ASSERT (cur_score <= n_col) ; - COLAMD_ASSERT (cur_score >= EIGEN_COLAMD_EMPTY) ; - if (prev_col == EIGEN_COLAMD_EMPTY) + COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ; + if (prev_col == COLAMD_EMPTY) { head [cur_score] = next_col ; } @@ -1379,7 +1225,7 @@ { Col [prev_col].shared4.degree_next = next_col ; } - if (next_col != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = prev_col ; } @@ -1394,7 +1240,7 @@ row = *cp++ ; row_mark = Row [row].shared2.mark ; /* skip if dead */ - if (EIGEN_ROW_IS_MARKED_DEAD (row_mark)) + if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } @@ -1413,7 +1259,7 @@ if (set_difference == 0) { COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; - EIGEN_KILL_ROW (row) ; + KILL_ROW (row) ; } else { @@ -1423,10 +1269,6 @@ } } -#ifndef COLAMD_NDEBUG - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k-pivot_row_degree, max_deg) ; -#endif /* COLAMD_NDEBUG */ /* === Add up set differences for each column ======================= */ @@ -1439,7 +1281,7 @@ { /* get a column */ col = *rp++ ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col) && col != pivot_col) ; + COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; hash = 0 ; cur_score = 0 ; cp = &A [Col [col].start] ; @@ -1456,7 +1298,7 @@ COLAMD_ASSERT(row >= 0 && row < n_row) ; row_mark = Row [row].shared2.mark ; /* skip if dead */ - if (EIGEN_ROW_IS_MARKED_DEAD (row_mark)) + if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } @@ -1480,7 +1322,7 @@ { COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ; /* nothing left but the pivot row in this column */ - EIGEN_KILL_PRINCIPAL_COL (col) ; + KILL_PRINCIPAL_COL (col) ; pivot_row_degree -= Col [col].shared1.thickness ; COLAMD_ASSERT (pivot_row_degree >= 0) ; /* order it */ @@ -1504,7 +1346,7 @@ COLAMD_ASSERT (hash <= n_col) ; head_column = head [hash] ; - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { /* degree list "hash" is non-empty, use prev (shared3) of */ /* first column in degree list as head of hash bucket */ @@ -1521,7 +1363,7 @@ /* save hash function in Col [col].shared3.hash */ Col [col].shared3.hash = (int) hash ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col)) ; + COLAMD_ASSERT (COL_IS_ALIVE (col)) ; } } @@ -1531,17 +1373,13 @@ COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ; - eigen_detect_super_cols ( - -#ifndef COLAMD_NDEBUG - n_col, Row, -#endif /* COLAMD_NDEBUG */ + detect_super_cols ( Col, A, head, pivot_row_start, pivot_row_length) ; /* === Kill the pivotal column ====================================== */ - EIGEN_KILL_PRINCIPAL_COL (pivot_col) ; + KILL_PRINCIPAL_COL (pivot_col) ; /* === Clear mark =================================================== */ @@ -1549,14 +1387,9 @@ if (tag_mark >= max_mark) { COLAMD_DEBUG2 (("clearing tag_mark\n")) ; - tag_mark = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; } -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG3 (("check3\n")) ; - eigen_debug_mark (n_row, Row, tag_mark, max_mark) ; -#endif /* COLAMD_NDEBUG */ - /* === Finalize the new pivot row, and column scores ================ */ COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ; @@ -1570,7 +1403,7 @@ { col = *rp++ ; /* skip dead columns */ - if (EIGEN_COL_IS_DEAD (col)) + if (COL_IS_DEAD (col)) { continue ; } @@ -1604,11 +1437,11 @@ COLAMD_ASSERT (min_score <= n_col) ; COLAMD_ASSERT (cur_score >= 0) ; COLAMD_ASSERT (cur_score <= n_col) ; - COLAMD_ASSERT (head [cur_score] >= EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ; next_col = head [cur_score] ; Col [col].shared4.degree_next = next_col ; - Col [col].shared3.prev = EIGEN_COLAMD_EMPTY ; - if (next_col != EIGEN_COLAMD_EMPTY) + Col [col].shared3.prev = COLAMD_EMPTY ; + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = col ; } @@ -1619,11 +1452,6 @@ } -#ifndef COLAMD_NDEBUG - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k, max_deg) ; -#endif /* COLAMD_NDEBUG */ - /* === Resurrect the new pivot row ================================== */ if (pivot_row_degree > 0) @@ -1645,11 +1473,11 @@ /* ========================================================================== */ -/* === eigen_order_children ======================================================= */ +/* === order_children ======================================================= */ /* ========================================================================== */ /* - The eigen_find_ordering routine has ordered all of the principal columns (the + 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 @@ -1661,12 +1489,12 @@ columns. Not user-callable. */ - void eigen_order_children +static inline void order_children ( /* === Parameters ======================================================= */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + colamd_col Col [], /* of size n_col+1 */ int p [] /* p [0 ... n_col-1] is the column permutation*/ ) { @@ -1682,15 +1510,15 @@ for (i = 0 ; i < n_col ; i++) { /* find an un-ordered non-principal column */ - COLAMD_ASSERT (EIGEN_COL_IS_DEAD (i)) ; - if (!EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EIGEN_COLAMD_EMPTY) + COLAMD_ASSERT (COL_IS_DEAD (i)) ; + if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY) { parent = i ; /* once found, find its principal parent */ do { parent = Col [parent].shared1.parent ; - } while (!EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL (parent)) ; + } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; /* now, order all un-ordered non-principal columns along path */ /* to this parent. collapse tree at the same time */ @@ -1700,7 +1528,7 @@ do { - COLAMD_ASSERT (Col [c].shared2.order == EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ; /* order this column */ Col [c].shared2.order = order++ ; @@ -1713,7 +1541,7 @@ /* continue until we hit an ordered column. There are */ /* guarranteed not to be anymore unordered columns */ /* above an ordered column */ - } while (Col [c].shared2.order == EIGEN_COLAMD_EMPTY) ; + } while (Col [c].shared2.order == COLAMD_EMPTY) ; /* re-order the super_col parent to largest order for this group */ Col [parent].shared2.order = order ; @@ -1730,7 +1558,7 @@ /* ========================================================================== */ -/* === eigen_detect_super_cols ==================================================== */ +/* === detect_super_cols ==================================================== */ /* ========================================================================== */ /* @@ -1762,17 +1590,11 @@ Not user-callable. */ - void eigen_detect_super_cols +static void detect_super_cols ( /* === Parameters ======================================================= */ -#ifndef COLAMD_NDEBUG - /* these two parameters are only needed when debugging is enabled: */ - int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ -#endif /* COLAMD_NDEBUG */ - - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* row indices of A */ int head [], /* head of degree lists and hash buckets */ int row_start, /* pointer to set of columns to check */ @@ -1802,7 +1624,7 @@ while (rp < rp_end) { col = *rp++ ; - if (EIGEN_COL_IS_DEAD (col)) + if (COL_IS_DEAD (col)) { continue ; } @@ -1814,7 +1636,7 @@ /* === Get the first column in this hash bucket ===================== */ head_column = head [hash] ; - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { first_col = Col [head_column].shared3.headhash ; } @@ -1825,10 +1647,10 @@ /* === Consider each column in the hash bucket ====================== */ - for (super_c = first_col ; super_c != EIGEN_COLAMD_EMPTY ; + for (super_c = first_col ; super_c != COLAMD_EMPTY ; super_c = Col [super_c].shared4.hash_next) { - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (super_c)) ; + COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ; COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ; length = Col [super_c].length ; @@ -1838,10 +1660,10 @@ /* === Compare super_c with all columns after it ================ */ for (c = Col [super_c].shared4.hash_next ; - c != EIGEN_COLAMD_EMPTY ; c = Col [c].shared4.hash_next) + c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next) { COLAMD_ASSERT (c != super_c) ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (c)) ; + COLAMD_ASSERT (COL_IS_ALIVE (c)) ; COLAMD_ASSERT (Col [c].shared3.hash == hash) ; /* not identical if lengths or scores are different */ @@ -1859,8 +1681,8 @@ for (i = 0 ; i < length ; i++) { /* the columns are "clean" (no dead rows) */ - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (*cp1)) ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (*cp2)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ; /* row indices will same order for both supercols, */ /* no gather scatter nessasary */ if (*cp1++ != *cp2++) @@ -1882,9 +1704,9 @@ Col [super_c].shared1.thickness += Col [c].shared1.thickness ; Col [c].shared1.parent = super_c ; - EIGEN_KILL_NON_PRINCIPAL_COL (c) ; - /* order c later, in eigen_order_children() */ - Col [c].shared2.order = EIGEN_COLAMD_EMPTY ; + KILL_NON_PRINCIPAL_COL (c) ; + /* order c later, in order_children() */ + Col [c].shared2.order = COLAMD_EMPTY ; /* remove c from hash bucket */ Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; } @@ -1892,22 +1714,22 @@ /* === Empty this hash bucket ======================================= */ - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { /* corresponding degree list "hash" is not empty */ - Col [head_column].shared3.headhash = EIGEN_COLAMD_EMPTY ; + Col [head_column].shared3.headhash = COLAMD_EMPTY ; } else { /* corresponding degree list "hash" is empty */ - head [hash] = EIGEN_COLAMD_EMPTY ; + head [hash] = COLAMD_EMPTY ; } } } /* ========================================================================== */ -/* === eigen_garbage_collection =================================================== */ +/* === garbage_collection =================================================== */ /* ========================================================================== */ /* @@ -1919,14 +1741,14 @@ Not user-callable. */ - int eigen_garbage_collection /* returns the new value of pfree */ +static int garbage_collection /* returns the new value of pfree */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows */ int n_col, /* number of columns */ - EIGEN_Colamd_Row Row [], /* row info */ - EIGEN_Colamd_Col Col [], /* column info */ + Colamd_Row Row [], /* row info */ + colamd_col Col [], /* column info */ int A [], /* A [0 ... Alen-1] holds the matrix */ int *pfree /* &A [0] ... pfree is in use */ ) @@ -1940,19 +1762,12 @@ int c ; /* a column index */ int length ; /* length of a row or column */ -#ifndef COLAMD_NDEBUG - int debug_rows ; - COLAMD_DEBUG2 (("Defrag..\n")) ; - for (psrc = &A[0] ; psrc < pfree ; psrc++) COLAMD_ASSERT (*psrc >= 0) ; - debug_rows = 0 ; -#endif /* COLAMD_NDEBUG */ - /* === Defragment the columns =========================================== */ pdest = &A[0] ; for (c = 0 ; c < n_col ; c++) { - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { psrc = &A [Col [c].start] ; @@ -1963,7 +1778,7 @@ for (j = 0 ; j < length ; j++) { r = *psrc++ ; - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { *pdest++ = r ; } @@ -1976,26 +1791,22 @@ for (r = 0 ; r < n_row ; r++) { - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { if (Row [r].length == 0) { /* this row is of zero length. cannot compact it, so kill it */ COLAMD_DEBUG3 (("Defrag row kill\n")) ; - EIGEN_KILL_ROW (r) ; + KILL_ROW (r) ; } else { /* save first column index in Row [r].shared2.first_column */ psrc = &A [Row [r].start] ; Row [r].shared2.first_column = *psrc ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; /* flag the start of the row with the one's complement of row */ - *psrc = EIGEN_ONES_COMPLEMENT (r) ; - -#ifndef COLAMD_NDEBUG - debug_rows++ ; -#endif /* COLAMD_NDEBUG */ + *psrc = ONES_COMPLEMENT (r) ; } } @@ -2011,11 +1822,11 @@ { psrc-- ; /* get the row index */ - r = EIGEN_ONES_COMPLEMENT (*psrc) ; + r = ONES_COMPLEMENT (*psrc) ; COLAMD_ASSERT (r >= 0 && r < n_row) ; /* restore first column index */ *psrc = Row [r].shared2.first_column ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; /* move and compact the row */ COLAMD_ASSERT (pdest <= psrc) ; @@ -2024,17 +1835,13 @@ for (j = 0 ; j < length ; j++) { c = *psrc++ ; - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { *pdest++ = c ; } } Row [r].length = (int) (pdest - &A [Row [r].start]) ; -#ifndef COLAMD_NDEBUG - debug_rows-- ; -#endif /* COLAMD_NDEBUG */ - } } /* ensure we found all the rows */ @@ -2047,7 +1854,7 @@ /* ========================================================================== */ -/* === eigen_clear_mark =========================================================== */ +/* === clear_mark =========================================================== */ /* ========================================================================== */ /* @@ -2055,12 +1862,12 @@ Return value is the new tag_mark. Not user-callable. */ - int eigen_clear_mark /* return the new value for tag_mark */ +static inline int clear_mark /* return the new value for tag_mark */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows in A */ - EIGEN_Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ + Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ ) { /* === Local variables ================================================== */ @@ -2069,7 +1876,7 @@ for (r = 0 ; r < n_row ; r++) { - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { Row [r].shared2.mark = 0 ; } @@ -2080,13 +1887,13 @@ /* ========================================================================== */ -/* === eigen_print_report ========================================================= */ +/* === print_report ========================================================= */ /* ========================================================================== */ - void eigen_print_report +static void print_report ( - char *method, - int stats [EIGEN_COLAMD_STATS] + const char *method, + int stats [COLAMD_STATS] ) { @@ -2098,11 +1905,11 @@ return ; } - i1 = stats [EIGEN_COLAMD_INFO1] ; - i2 = stats [EIGEN_COLAMD_INFO2] ; - i3 = stats [EIGEN_COLAMD_INFO3] ; + i1 = stats [COLAMD_INFO1] ; + i2 = stats [COLAMD_INFO2] ; + i3 = stats [COLAMD_INFO3] ; - if (stats [EIGEN_COLAMD_STATUS] >= 0) + if (stats [COLAMD_STATUS] >= 0) { PRINTF ("%s: OK. ", method) ; } @@ -2111,10 +1918,10 @@ PRINTF ("%s: ERROR. ", method) ; } - switch (stats [EIGEN_COLAMD_STATUS]) + switch (stats [COLAMD_STATUS]) { - case EIGEN_COLAMD_OK_BUT_JUMBLED: + case COLAMD_OK_BUT_JUMBLED: PRINTF ("Matrix has unsorted or duplicate row indices.\n") ; @@ -2129,77 +1936,77 @@ /* no break - fall through to next case instead */ - case EIGEN_COLAMD_OK: + case COLAMD_OK: PRINTF ("\n") ; PRINTF ("%s: number of dense or empty rows ignored: %d\n", - method, stats [EIGEN_COLAMD_DENSE_ROW]) ; + method, stats [COLAMD_DENSE_ROW]) ; PRINTF ("%s: number of dense or empty columns ignored: %d\n", - method, stats [EIGEN_COLAMD_DENSE_COL]) ; + method, stats [COLAMD_DENSE_COL]) ; PRINTF ("%s: number of garbage collections performed: %d\n", - method, stats [EIGEN_COLAMD_DEFRAG_COUNT]) ; + method, stats [COLAMD_DEFRAG_COUNT]) ; break ; - case EIGEN_COLAMD_ERROR_A_not_present: + case COLAMD_ERROR_A_not_present: PRINTF ("Array A (row indices of matrix) not present.\n") ; break ; - case EIGEN_COLAMD_ERROR_p_not_present: + case COLAMD_ERROR_p_not_present: PRINTF ("Array p (column pointers for matrix) not present.\n") ; break ; - case EIGEN_COLAMD_ERROR_nrow_negative: + case COLAMD_ERROR_nrow_negative: PRINTF ("Invalid number of rows (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_ncol_negative: + case COLAMD_ERROR_ncol_negative: PRINTF ("Invalid number of columns (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_nnz_negative: + case COLAMD_ERROR_nnz_negative: PRINTF ("Invalid number of nonzero entries (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_p0_nonzero: + case COLAMD_ERROR_p0_nonzero: PRINTF ("Invalid column pointer, p [0] = %d, must be zero.\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_A_too_small: + case COLAMD_ERROR_A_too_small: PRINTF ("Array A too small.\n") ; PRINTF (" Need Alen >= %d, but given only Alen = %d.\n", i1, i2) ; break ; - case EIGEN_COLAMD_ERROR_col_length_negative: + case COLAMD_ERROR_col_length_negative: PRINTF ("Column %d has a negative number of nonzero entries (%d).\n", INDEX (i1), i2) ; break ; - case EIGEN_COLAMD_ERROR_row_index_out_of_bounds: + case COLAMD_ERROR_row_index_out_of_bounds: PRINTF ("Row index (row %d) out of bounds (%d to %d) in column %d.\n", INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1)) ; break ; - case EIGEN_COLAMD_ERROR_out_of_memory: + case COLAMD_ERROR_out_of_memory: PRINTF ("Out of memory.\n") ; break ; - case EIGEN_COLAMD_ERROR_internal_error: + case COLAMD_ERROR_internal_error: /* if this happens, there is a bug in the code */ PRINTF @@ -2208,307 +2015,5 @@ } } - - - -/* ========================================================================== */ -/* === eigen_colamd debugging routines ============================================ */ -/* ========================================================================== */ - -/* When debugging is disabled, the remainder of this file is ignored. */ - -#ifndef COLAMD_NDEBUG - - -/* ========================================================================== */ -/* === eigen_debug_structures ===================================================== */ -/* ========================================================================== */ - -/* - 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. -*/ - - void eigen_debug_structures -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [], - int n_col2 -) -{ - /* === Local variables ================================================== */ - - int i ; - int c ; - int *cp ; - int *cp_end ; - int len ; - int score ; - int r ; - int *rp ; - int *rp_end ; - int deg ; - - /* === Check A, Row, and Col ============================================ */ - - for (c = 0 ; c < n_col ; c++) - { - if (EIGEN_COL_IS_ALIVE (c)) - { - len = Col [c].length ; - score = Col [c].shared2.score ; - COLAMD_DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; - COLAMD_ASSERT (len > 0) ; - COLAMD_ASSERT (score >= 0) ; - COLAMD_ASSERT (Col [c].shared1.thickness == 1) ; - cp = &A [Col [c].start] ; - cp_end = cp + len ; - while (cp < cp_end) - { - r = *cp++ ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; - } - } - else - { - i = Col [c].shared2.order ; - COLAMD_ASSERT (i >= n_col2 && i < n_col) ; - } - } - - for (r = 0 ; r < n_row ; r++) - { - if (EIGEN_ROW_IS_ALIVE (r)) - { - i = 0 ; - len = Row [r].length ; - deg = Row [r].shared1.degree ; - COLAMD_ASSERT (len > 0) ; - COLAMD_ASSERT (deg > 0) ; - rp = &A [Row [r].start] ; - rp_end = rp + len ; - while (rp < rp_end) - { - c = *rp++ ; - if (EIGEN_COL_IS_ALIVE (c)) - { - i++ ; - } - } - COLAMD_ASSERT (i > 0) ; - } - } -} - - -/* ========================================================================== */ -/* === eigen_debug_deg_lists ====================================================== */ -/* ========================================================================== */ - -/* - Prints the contents of the degree lists. Counts the number of columns - in the degree list and compares it to the total it should have. Also - checks the row degrees. -*/ - - void eigen_debug_deg_lists -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int head [], - int min_score, - int should, - int max_deg -) -{ - /* === Local variables ================================================== */ - - int deg ; - int col ; - int have ; - int row ; - - /* === Check the degree lists =========================================== */ - - if (n_col > 10000 && colamd_debug <= 0) - { - return ; - } - have = 0 ; - COLAMD_DEBUG4 (("Degree lists: %d\n", min_score)) ; - for (deg = 0 ; deg <= n_col ; deg++) - { - col = head [deg] ; - if (col == EIGEN_COLAMD_EMPTY) - { - continue ; - } - COLAMD_DEBUG4 (("%d:", deg)) ; - while (col != EIGEN_COLAMD_EMPTY) - { - COLAMD_DEBUG4 ((" %d", col)) ; - have += Col [col].shared1.thickness ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col)) ; - col = Col [col].shared4.degree_next ; - } - COLAMD_DEBUG4 (("\n")) ; - } - COLAMD_DEBUG4 (("should %d have %d\n", should, have)) ; - COLAMD_ASSERT (should == have) ; - - /* === Check the row degrees ============================================ */ - - if (n_row > 10000 && colamd_debug <= 0) - { - return ; - } - for (row = 0 ; row < n_row ; row++) - { - if (EIGEN_ROW_IS_ALIVE (row)) - { - COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ; - } - } -} - - -/* ========================================================================== */ -/* === eigen_debug_mark =========================================================== */ -/* ========================================================================== */ - -/* - Ensures that the tag_mark is less that the maximum and also ensures that - each entry in the mark array is less than the tag mark. -*/ - - void eigen_debug_mark -( - /* === Parameters ======================================================= */ - - int n_row, - EIGEN_Colamd_Row Row [], - int tag_mark, - int max_mark -) -{ - /* === Local variables ================================================== */ - - int r ; - - /* === Check the Row marks ============================================== */ - - COLAMD_ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; - if (n_row > 10000 && colamd_debug <= 0) - { - return ; - } - for (r = 0 ; r < n_row ; r++) - { - COLAMD_ASSERT (Row [r].shared2.mark < tag_mark) ; - } -} - - -/* ========================================================================== */ -/* === eigen_debug_matrix ========================================================= */ -/* ========================================================================== */ - -/* - Prints out the contents of the columns and the rows. -*/ - - void eigen_debug_matrix -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [] -) -{ - /* === Local variables ================================================== */ - - int r ; - int c ; - int *rp ; - int *rp_end ; - int *cp ; - int *cp_end ; - - /* === Dump the rows and columns of the matrix ========================== */ - - if (colamd_debug < 3) - { - return ; - } - COLAMD_DEBUG3 (("DUMP MATRIX:\n")) ; - for (r = 0 ; r < n_row ; r++) - { - COLAMD_DEBUG3 (("Row %d alive? %d\n", r, EIGEN_ROW_IS_ALIVE (r))) ; - if (EIGEN_ROW_IS_DEAD (r)) - { - continue ; - } - COLAMD_DEBUG3 (("start %d length %d degree %d\n", - Row [r].start, Row [r].length, Row [r].shared1.degree)) ; - rp = &A [Row [r].start] ; - rp_end = rp + Row [r].length ; - while (rp < rp_end) - { - c = *rp++ ; - COLAMD_DEBUG4 ((" %d col %d\n", EIGEN_COL_IS_ALIVE (c), c)) ; - } - } - - for (c = 0 ; c < n_col ; c++) - { - COLAMD_DEBUG3 (("Col %d alive? %d\n", c, EIGEN_COL_IS_ALIVE (c))) ; - if (EIGEN_COL_IS_DEAD (c)) - { - continue ; - } - COLAMD_DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", - Col [c].start, Col [c].length, - Col [c].shared1.thickness, Col [c].shared2.score)) ; - cp = &A [Col [c].start] ; - cp_end = cp + Col [c].length ; - while (cp < cp_end) - { - r = *cp++ ; - COLAMD_DEBUG4 ((" %d row %d\n", EIGEN_ROW_IS_ALIVE (r), r)) ; - } - } -} - - void eigen_colamd_get_debug -( - char *method -) -{ - colamd_debug = 0 ; /* no debug printing */ - - /* get "D" environment variable, which gives the debug printing level */ - if (getenv ("D")) - { - colamd_debug = atoi (getenv ("D")) ; - } - - COLAMD_DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", - method, colamd_debug)) ; -} - -#endif /* NDEBUG */ +} // namespace internal #endif
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index 47cd6f1..f5757b3 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -27,8 +27,10 @@ #define EIGEN_ORDERING_H #include "Amd.h" -#include "Eigen_Colamd.h" namespace Eigen { + +#include "Eigen_Colamd.h" + namespace internal { /** @@ -131,18 +133,18 @@ int n = mat.cols(); int nnz = mat.nonZeros(); // Get the recommended value of Alen to be used by colamd - int Alen = eigen_colamd_recommended(nnz, m, n); + int Alen = internal::colamd_recommended(nnz, m, n); // Set the default parameters - double knobs [EIGEN_COLAMD_KNOBS]; - int stats [EIGEN_COLAMD_STATS]; - eigen_colamd_set_defaults(knobs); + double knobs [COLAMD_KNOBS]; + int stats [COLAMD_STATS]; + internal::colamd_set_defaults(knobs); int info; IndexVector p(n+1), A(Alen); for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; // Call Colamd routine to compute the ordering - info = eigen_colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); eigen_assert( info && "COLAMD failed " ); perm.resize(n);
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index e207613..77df091 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h
@@ -205,7 +205,7 @@ void initperfvalues() { m_perfv.panel_size = 12; - m_perfv.relax = 6; + m_perfv.relax = 1; m_perfv.maxsuper = 100; m_perfv.rowblk = 200; m_perfv.colblk = 60;
diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 5451843..6e0e1b1 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt
@@ -55,6 +55,12 @@ set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES} ${BLAS_LIBRARIES}) endif(PASTIX_FOUND AND BLAS_FOUND) +if(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) + add_definitions("-DEIGEN_METIS_SUPPORT") +endif(METIS_FOUND) + find_library(RT_LIBRARY rt) if(RT_LIBRARY) set(SPARSE_LIBS ${SPARSE_LIBS} ${RT_LIBRARY}) @@ -66,11 +72,6 @@ add_executable(spsolver sp_solver.cpp) target_link_libraries (spsolver ${SPARSE_LIBS}) -if(METIS_FOUND) - include_directories(${METIS_INCLUDES}) - set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) - add_definitions("-DEIGEN_METIS_SUPPORT") -endif(METIS_FOUND) add_executable(test_sparseLU test_sparseLU.cpp) target_link_libraries (test_sparseLU ${SPARSE_LIBS})
diff --git a/bench/spbench/spbenchsolver.h b/bench/spbench/spbenchsolver.h index c48ed7a..19c719c 100644 --- a/bench/spbench/spbenchsolver.h +++ b/bench/spbench/spbenchsolver.h
@@ -21,9 +21,14 @@ #include <unsupported/Eigen/IterativeSolvers> #include <Eigen/LU> #include <unsupported/Eigen/SparseExtra> +#include <Eigen/SparseLU> #include "spbenchstyle.h" +#ifdef EIGEN_METIS_SUPPORT +#include <Eigen/MetisSupport> +#endif + #ifdef EIGEN_CHOLMOD_SUPPORT #include <Eigen/CholmodSupport> #endif @@ -45,26 +50,27 @@ #endif // CONSTANTS -#define EIGEN_UMFPACK 0 -#define EIGEN_SUPERLU 1 -#define EIGEN_PASTIX 2 -#define EIGEN_PARDISO 3 -#define EIGEN_BICGSTAB 4 -#define EIGEN_BICGSTAB_ILUT 5 -#define EIGEN_GMRES 6 -#define EIGEN_GMRES_ILUT 7 -#define EIGEN_SIMPLICIAL_LDLT 8 -#define EIGEN_CHOLMOD_LDLT 9 -#define EIGEN_PASTIX_LDLT 10 -#define EIGEN_PARDISO_LDLT 11 -#define EIGEN_SIMPLICIAL_LLT 12 -#define EIGEN_CHOLMOD_SUPERNODAL_LLT 13 -#define EIGEN_CHOLMOD_SIMPLICIAL_LLT 14 -#define EIGEN_PASTIX_LLT 15 -#define EIGEN_PARDISO_LLT 16 -#define EIGEN_CG 17 -#define EIGEN_CG_PRECOND 18 -#define EIGEN_ALL_SOLVERS 19 +#define EIGEN_UMFPACK 10 +#define EIGEN_SUPERLU 20 +#define EIGEN_PASTIX 30 +#define EIGEN_PARDISO 40 +#define EIGEN_SPARSELU_COLAMD 50 +#define EIGEN_SPARSELU_METIS 51 +#define EIGEN_BICGSTAB 60 +#define EIGEN_BICGSTAB_ILUT 61 +#define EIGEN_GMRES 70 +#define EIGEN_GMRES_ILUT 71 +#define EIGEN_SIMPLICIAL_LDLT 80 +#define EIGEN_CHOLMOD_LDLT 90 +#define EIGEN_PASTIX_LDLT 100 +#define EIGEN_PARDISO_LDLT 110 +#define EIGEN_SIMPLICIAL_LLT 120 +#define EIGEN_CHOLMOD_SUPERNODAL_LLT 130 +#define EIGEN_CHOLMOD_SIMPLICIAL_LLT 140 +#define EIGEN_PASTIX_LLT 150 +#define EIGEN_PARDISO_LLT 160 +#define EIGEN_CG 170 +#define EIGEN_CG_PRECOND 180 using namespace Eigen; using namespace std; @@ -188,6 +194,17 @@ out << " <PACKAGE> EIGEN </PACKAGE> \n"; out << " </SOLVER> \n"; + out <<" <SOLVER ID='" << EIGEN_SPARSELU_COLAMD << "'>\n"; + out << " <TYPE> LU_COLAMD </TYPE> \n"; + out << " <PACKAGE> EIGEN </PACKAGE> \n"; + out << " </SOLVER> \n"; + +#ifdef EIGEN_METIS_SUPPORT + out <<" <SOLVER ID='" << EIGEN_SPARSELU_METIS << "'>\n"; + out << " <TYPE> LU_METIS </TYPE> \n"; + out << " <PACKAGE> EIGEN </PACKAGE> \n"; + out << " </SOLVER> \n"; +#endif out << " </AVAILSOLVER> \n"; } @@ -325,8 +342,19 @@ call_directsolver(solver, EIGEN_PARDISO, A, b, refX,statFile); } #endif - - + + // Eigen SparseLU METIS + cout << "\n Solving with Sparse LU AND COLAMD ... \n"; + SparseLU<SpMat, COLAMDOrdering<int> > solver; + call_directsolver(solver, EIGEN_SPARSELU_COLAMD, A, b, refX, statFile); + // Eigen SparseLU METIS + #ifdef EIGEN_METIS_SUPPORT + { + cout << "\n Solving with Sparse LU AND METIS ... \n"; + SparseLU<SpMat, MetisOrdering<int> > solver; + call_directsolver(solver, EIGEN_SPARSELU_METIS, A, b, refX, statFile); + } + #endif //BiCGSTAB {