47 #ifndef EIGEN_COLAMD_H 48 #define EIGEN_COLAMD_H 60 #define COLAMD_KNOBS 20 63 #define COLAMD_STATS 20 66 #define COLAMD_DENSE_ROW 0 69 #define COLAMD_DENSE_COL 1 72 #define COLAMD_DEFRAG_COUNT 2 75 #define COLAMD_STATUS 3 78 #define COLAMD_INFO1 4 79 #define COLAMD_INFO2 5 80 #define COLAMD_INFO3 6 84 #define COLAMD_OK_BUT_JUMBLED (1) 85 #define COLAMD_ERROR_A_not_present (-1) 86 #define COLAMD_ERROR_p_not_present (-2) 87 #define COLAMD_ERROR_nrow_negative (-3) 88 #define COLAMD_ERROR_ncol_negative (-4) 89 #define COLAMD_ERROR_nnz_negative (-5) 90 #define COLAMD_ERROR_p0_nonzero (-6) 91 #define COLAMD_ERROR_A_too_small (-7) 92 #define COLAMD_ERROR_col_length_negative (-8) 93 #define COLAMD_ERROR_row_index_out_of_bounds (-9) 94 #define COLAMD_ERROR_out_of_memory (-10) 95 #define COLAMD_ERROR_internal_error (-999) 101 #define ONES_COMPLEMENT(r) (-(r)-1) 105 #define COLAMD_EMPTY (-1) 112 #define DEAD_PRINCIPAL (-1) 113 #define DEAD_NON_PRINCIPAL (-2) 116 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) 117 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) 118 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) 119 #define COL_IS_DEAD(c) (Col [c].start < ALIVE) 120 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) 121 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) 122 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } 123 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } 124 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } 131 template <
typename IndexType>
139 IndexType thickness ;
159 IndexType degree_next ;
160 IndexType hash_next ;
165 template <
typename IndexType>
178 IndexType first_column ;
201 template <
typename IndexType>
202 inline IndexType colamd_c(IndexType n_col)
203 {
return IndexType( ((n_col) + 1) *
sizeof (colamd_col<IndexType>) /
sizeof (IndexType) ) ; }
205 template <
typename IndexType>
206 inline IndexType colamd_r(IndexType n_row)
207 {
return IndexType(((n_row) + 1) *
sizeof (Colamd_Row<IndexType>) /
sizeof (IndexType)); }
210 template <
typename IndexType>
211 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] );
213 template <
typename IndexType>
214 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);
216 template <
typename IndexType>
217 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);
219 template <
typename IndexType>
220 static void order_children (IndexType n_col, colamd_col<IndexType> Col [], IndexType p []);
222 template <
typename IndexType>
223 static void detect_super_cols (colamd_col<IndexType> Col [], IndexType A [], IndexType head [], IndexType row_start, IndexType row_length ) ;
225 template <
typename IndexType>
226 static IndexType garbage_collection (IndexType n_row, IndexType n_col, Colamd_Row<IndexType> Row [], colamd_col<IndexType> Col [], IndexType A [], IndexType *pfree) ;
228 template <
typename IndexType>
229 static inline IndexType clear_mark (IndexType n_row, Colamd_Row<IndexType> Row [] ) ;
233 #define COLAMD_DEBUG0(params) ; 234 #define COLAMD_DEBUG1(params) ; 235 #define COLAMD_DEBUG2(params) ; 236 #define COLAMD_DEBUG3(params) ; 237 #define COLAMD_DEBUG4(params) ; 239 #define COLAMD_ASSERT(expression) ((void) 0) 256 template <
typename IndexType>
257 inline IndexType colamd_recommended ( IndexType nnz, IndexType n_row, IndexType n_col)
259 if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0)
262 return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5));
286 static inline void colamd_set_defaults(
double knobs[COLAMD_KNOBS])
296 for (i = 0 ; i < COLAMD_KNOBS ; i++)
300 knobs [COLAMD_DENSE_ROW] = 0.5 ;
301 knobs [COLAMD_DENSE_COL] = 0.5 ;
321 template <
typename IndexType>
322 static bool colamd(IndexType n_row, IndexType n_col, IndexType Alen, IndexType *A, IndexType *p,
double knobs[COLAMD_KNOBS], IndexType stats[COLAMD_STATS])
331 Colamd_Row<IndexType> *Row ;
332 colamd_col<IndexType> *Col ;
337 double default_knobs [COLAMD_KNOBS] ;
344 COLAMD_DEBUG0 ((
"colamd: stats not present\n")) ;
347 for (i = 0 ; i < COLAMD_STATS ; i++)
351 stats [COLAMD_STATUS] = COLAMD_OK ;
352 stats [COLAMD_INFO1] = -1 ;
353 stats [COLAMD_INFO2] = -1 ;
357 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ;
358 COLAMD_DEBUG0 ((
"colamd: A not present\n")) ;
364 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ;
365 COLAMD_DEBUG0 ((
"colamd: p not present\n")) ;
371 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ;
372 stats [COLAMD_INFO1] = n_row ;
373 COLAMD_DEBUG0 ((
"colamd: nrow negative %d\n", n_row)) ;
379 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ;
380 stats [COLAMD_INFO1] = n_col ;
381 COLAMD_DEBUG0 ((
"colamd: ncol negative %d\n", n_col)) ;
388 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ;
389 stats [COLAMD_INFO1] = nnz ;
390 COLAMD_DEBUG0 ((
"colamd: number of entries negative %d\n", nnz)) ;
396 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ;
397 stats [COLAMD_INFO1] = p [0] ;
398 COLAMD_DEBUG0 ((
"colamd: p[0] not zero %d\n", p [0])) ;
406 colamd_set_defaults (default_knobs) ;
407 knobs = default_knobs ;
412 Col_size = colamd_c (n_col) ;
413 Row_size = colamd_r (n_row) ;
414 need = 2*nnz + n_col + Col_size + Row_size ;
419 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ;
420 stats [COLAMD_INFO1] = need ;
421 stats [COLAMD_INFO2] = Alen ;
422 COLAMD_DEBUG0 ((
"colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
426 Alen -= Col_size + Row_size ;
427 Col = (colamd_col<IndexType> *) &A [Alen] ;
428 Row = (Colamd_Row<IndexType> *) &A [Alen + Col_size] ;
432 if (!Eigen::internal::init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
435 COLAMD_DEBUG0 ((
"colamd: Matrix invalid\n")) ;
441 Eigen::internal::init_scoring (n_row, n_col, Row, Col, A, p, knobs,
442 &n_row2, &n_col2, &max_deg) ;
446 ngarbage = Eigen::internal::find_ordering (n_row, n_col, Alen, Row, Col, A, p,
447 n_col2, max_deg, 2*nnz) ;
451 Eigen::internal::order_children (n_col, Col, p) ;
455 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ;
456 stats [COLAMD_DENSE_COL] = n_col - n_col2 ;
457 stats [COLAMD_DEFRAG_COUNT] = ngarbage ;
458 COLAMD_DEBUG0 ((
"colamd: done.\n")) ;
481 template <
typename IndexType>
482 static IndexType init_rows_cols
488 Colamd_Row<IndexType> Row [],
489 colamd_col<IndexType> Col [],
492 IndexType stats [COLAMD_STATS]
507 for (col = 0 ; col < n_col ; col++)
509 Col [col].start = p [col] ;
510 Col [col].length = p [col+1] - p [col] ;
512 if ((Col [col].length) < 0)
515 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ;
516 stats [COLAMD_INFO1] = col ;
517 stats [COLAMD_INFO2] = Col [col].length ;
518 COLAMD_DEBUG0 ((
"colamd: col %d length %d < 0\n", col, Col [col].length)) ;
522 Col [col].shared1.thickness = 1 ;
523 Col [col].shared2.score = 0 ;
524 Col [col].shared3.prev = COLAMD_EMPTY ;
525 Col [col].shared4.degree_next = COLAMD_EMPTY ;
532 stats [COLAMD_INFO3] = 0 ;
534 for (row = 0 ; row < n_row ; row++)
536 Row [row].length = 0 ;
537 Row [row].shared2.mark = -1 ;
540 for (col = 0 ; col < n_col ; col++)
545 cp_end = &A [p [col+1]] ;
552 if (row < 0 || row >= n_row)
554 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ;
555 stats [COLAMD_INFO1] = col ;
556 stats [COLAMD_INFO2] = row ;
557 stats [COLAMD_INFO3] = n_row ;
558 COLAMD_DEBUG0 ((
"colamd: row %d col %d out of bounds\n", row, col)) ;
562 if (row <= last_row || Row [row].shared2.mark == col)
566 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ;
567 stats [COLAMD_INFO1] = col ;
568 stats [COLAMD_INFO2] = row ;
569 (stats [COLAMD_INFO3]) ++ ;
570 COLAMD_DEBUG1 ((
"colamd: row %d col %d unsorted/duplicate\n",row,col));
573 if (Row [row].shared2.mark != col)
585 Row [row].shared2.mark = col ;
595 Row [0].start = p [n_col] ;
596 Row [0].shared1.p = Row [0].start ;
597 Row [0].shared2.mark = -1 ;
598 for (row = 1 ; row < n_row ; row++)
600 Row [row].start = Row [row-1].start + Row [row-1].length ;
601 Row [row].shared1.p = Row [row].start ;
602 Row [row].shared2.mark = -1 ;
607 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
610 for (col = 0 ; col < n_col ; col++)
613 cp_end = &A [p [col+1]] ;
617 if (Row [row].shared2.mark != col)
619 A [(Row [row].shared1.p)++] = col ;
620 Row [row].shared2.mark = col ;
628 for (col = 0 ; col < n_col ; col++)
631 cp_end = &A [p [col+1]] ;
634 A [(Row [*cp++].shared1.p)++] = col ;
641 for (row = 0 ; row < n_row ; row++)
643 Row [row].shared2.mark = 0 ;
644 Row [row].shared1.degree = Row [row].length ;
649 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED)
651 COLAMD_DEBUG0 ((
"colamd: reconstructing column form, matrix jumbled\n")) ;
661 p [0] = Col [0].start ;
662 for (col = 1 ; col < n_col ; col++)
666 Col [col].start = Col [col-1].start + Col [col-1].length ;
667 p [col] = Col [col].start ;
672 for (row = 0 ; row < n_row ; row++)
674 rp = &A [Row [row].start] ;
675 rp_end = rp + Row [row].length ;
678 A [(p [*rp++])++] = row ;
697 template <
typename IndexType>
698 static void init_scoring
704 Colamd_Row<IndexType> Row [],
705 colamd_col<IndexType> Col [],
708 double knobs [COLAMD_KNOBS],
722 IndexType col_length ;
726 IndexType dense_row_count ;
727 IndexType dense_col_count ;
728 IndexType min_score ;
735 dense_row_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_ROW] * n_col), n_col)) ;
736 dense_col_count = numext::maxi(IndexType(0), numext::mini(IndexType(knobs [COLAMD_DENSE_COL] * n_row), n_row)) ;
737 COLAMD_DEBUG1 ((
"colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
746 for (c = n_col-1 ; c >= 0 ; c--)
748 deg = Col [c].length ;
752 Col [c].shared2.order = --n_col2 ;
753 KILL_PRINCIPAL_COL (c) ;
756 COLAMD_DEBUG1 ((
"colamd: null columns killed: %d\n", n_col - n_col2)) ;
761 for (c = n_col-1 ; c >= 0 ; c--)
768 deg = Col [c].length ;
769 if (deg > dense_col_count)
772 Col [c].shared2.order = --n_col2 ;
774 cp = &A [Col [c].start] ;
775 cp_end = cp + Col [c].length ;
778 Row [*cp++].shared1.degree-- ;
780 KILL_PRINCIPAL_COL (c) ;
783 COLAMD_DEBUG1 ((
"colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
787 for (r = 0 ; r < n_row ; r++)
789 deg = Row [r].shared1.degree ;
790 COLAMD_ASSERT (deg >= 0 && deg <= n_col) ;
791 if (deg > dense_row_count || deg == 0)
800 max_deg = numext::maxi(max_deg, deg) ;
803 COLAMD_DEBUG1 ((
"colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
813 for (c = n_col-1 ; c >= 0 ; c--)
821 cp = &A [Col [c].start] ;
823 cp_end = cp + Col [c].length ;
829 if (ROW_IS_DEAD (row))
836 score += Row [row].shared1.degree - 1 ;
838 score = numext::mini(score, n_col) ;
841 col_length = (IndexType) (new_cp - &A [Col [c].start]) ;
846 COLAMD_DEBUG2 ((
"Newly null killed: %d\n", c)) ;
847 Col [c].shared2.order = --n_col2 ;
848 KILL_PRINCIPAL_COL (c) ;
853 COLAMD_ASSERT (score >= 0) ;
854 COLAMD_ASSERT (score <= n_col) ;
855 Col [c].length = col_length ;
856 Col [c].shared2.score = score ;
859 COLAMD_DEBUG1 ((
"colamd: Dense, null, and newly-null columns killed: %d\n",
871 for (c = 0 ; c <= n_col ; c++)
873 head [c] = COLAMD_EMPTY ;
878 for (c = n_col-1 ; c >= 0 ; c--)
881 if (COL_IS_ALIVE (c))
883 COLAMD_DEBUG4 ((
"place %d score %d minscore %d ncol %d\n",
884 c, Col [c].shared2.score, min_score, n_col)) ;
888 score = Col [c].shared2.score ;
890 COLAMD_ASSERT (min_score >= 0) ;
891 COLAMD_ASSERT (min_score <= n_col) ;
892 COLAMD_ASSERT (score >= 0) ;
893 COLAMD_ASSERT (score <= n_col) ;
894 COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ;
897 next_col = head [score] ;
898 Col [c].shared3.prev = COLAMD_EMPTY ;
899 Col [c].shared4.degree_next = next_col ;
903 if (next_col != COLAMD_EMPTY)
905 Col [next_col].shared3.prev = c ;
910 min_score = numext::mini(min_score, score) ;
921 *p_max_deg = max_deg ;
934 template <
typename IndexType>
935 static IndexType find_ordering
942 Colamd_Row<IndexType> Row [],
943 colamd_col<IndexType> Col [],
954 IndexType pivot_col ;
957 IndexType pivot_row ;
960 IndexType pivot_row_start ;
961 IndexType pivot_row_degree ;
962 IndexType pivot_row_length ;
963 IndexType pivot_col_score ;
964 IndexType needed_memory ;
969 IndexType max_score ;
970 IndexType cur_score ;
972 IndexType head_column ;
973 IndexType first_col ;
976 IndexType set_difference ;
977 IndexType min_score ;
978 IndexType col_thickness ;
980 IndexType pivot_col_thickness ;
988 max_mark = INT_MAX - n_col ;
989 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
992 COLAMD_DEBUG1 ((
"colamd: Ordering, n_col2=%d\n", n_col2)) ;
996 for (k = 0 ; k < n_col2 ; )
1002 COLAMD_ASSERT (min_score >= 0) ;
1003 COLAMD_ASSERT (min_score <= n_col) ;
1004 COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ;
1007 while (min_score < n_col && head [min_score] == COLAMD_EMPTY)
1011 pivot_col = head [min_score] ;
1012 COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
1013 next_col = Col [pivot_col].shared4.degree_next ;
1014 head [min_score] = next_col ;
1015 if (next_col != COLAMD_EMPTY)
1017 Col [next_col].shared3.prev = COLAMD_EMPTY ;
1020 COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ;
1021 COLAMD_DEBUG3 ((
"Pivot col: %d\n", pivot_col)) ;
1024 pivot_col_score = Col [pivot_col].shared2.score ;
1027 Col [pivot_col].shared2.order = k ;
1030 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
1031 k += pivot_col_thickness ;
1032 COLAMD_ASSERT (pivot_col_thickness > 0) ;
1036 needed_memory = numext::mini(pivot_col_score, n_col - k) ;
1037 if (pfree + needed_memory >= Alen)
1039 pfree = Eigen::internal::garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
1042 COLAMD_ASSERT (pfree + needed_memory < Alen) ;
1044 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
1051 pivot_row_start = pfree ;
1054 pivot_row_degree = 0 ;
1058 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
1061 cp = &A [Col [pivot_col].start] ;
1062 cp_end = cp + Col [pivot_col].length ;
1067 COLAMD_DEBUG4 ((
"Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ;
1069 if (ROW_IS_DEAD (row))
1073 rp = &A [Row [row].start] ;
1074 rp_end = rp + Row [row].length ;
1080 col_thickness = Col [col].shared1.thickness ;
1081 if (col_thickness > 0 && COL_IS_ALIVE (col))
1084 Col [col].shared1.thickness = -col_thickness ;
1085 COLAMD_ASSERT (pfree < Alen) ;
1088 pivot_row_degree += col_thickness ;
1094 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
1095 max_deg = numext::maxi(max_deg, pivot_row_degree) ;
1101 cp = &A [Col [pivot_col].start] ;
1102 cp_end = cp + Col [pivot_col].length ;
1107 COLAMD_DEBUG3 ((
"Kill row in pivot col: %d\n", row)) ;
1113 pivot_row_length = pfree - pivot_row_start ;
1114 if (pivot_row_length > 0)
1117 pivot_row = A [Col [pivot_col].start] ;
1118 COLAMD_DEBUG3 ((
"Pivotal row is %d\n", pivot_row)) ;
1123 pivot_row = COLAMD_EMPTY ;
1124 COLAMD_ASSERT (pivot_row_length == 0) ;
1126 COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
1149 COLAMD_DEBUG3 ((
"** Computing set differences phase. **\n")) ;
1153 COLAMD_DEBUG3 ((
"Pivot row: ")) ;
1155 rp = &A [pivot_row_start] ;
1156 rp_end = rp + pivot_row_length ;
1160 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1161 COLAMD_DEBUG3 ((
"Col: %d\n", col)) ;
1164 col_thickness = -Col [col].shared1.thickness ;
1165 COLAMD_ASSERT (col_thickness > 0) ;
1166 Col [col].shared1.thickness = col_thickness ;
1170 cur_score = Col [col].shared2.score ;
1171 prev_col = Col [col].shared3.prev ;
1172 next_col = Col [col].shared4.degree_next ;
1173 COLAMD_ASSERT (cur_score >= 0) ;
1174 COLAMD_ASSERT (cur_score <= n_col) ;
1175 COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ;
1176 if (prev_col == COLAMD_EMPTY)
1178 head [cur_score] = next_col ;
1182 Col [prev_col].shared4.degree_next = next_col ;
1184 if (next_col != COLAMD_EMPTY)
1186 Col [next_col].shared3.prev = prev_col ;
1191 cp = &A [Col [col].start] ;
1192 cp_end = cp + Col [col].length ;
1197 row_mark = Row [row].shared2.mark ;
1199 if (ROW_IS_MARKED_DEAD (row_mark))
1203 COLAMD_ASSERT (row != pivot_row) ;
1204 set_difference = row_mark - tag_mark ;
1206 if (set_difference < 0)
1208 COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ;
1209 set_difference = Row [row].shared1.degree ;
1212 set_difference -= col_thickness ;
1213 COLAMD_ASSERT (set_difference >= 0) ;
1215 if (set_difference == 0)
1217 COLAMD_DEBUG3 ((
"aggressive absorption. Row: %d\n", row)) ;
1223 Row [row].shared2.mark = set_difference + tag_mark ;
1231 COLAMD_DEBUG3 ((
"** Adding set differences phase. **\n")) ;
1234 rp = &A [pivot_row_start] ;
1235 rp_end = rp + pivot_row_length ;
1240 COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
1243 cp = &A [Col [col].start] ;
1246 cp_end = cp + Col [col].length ;
1248 COLAMD_DEBUG4 ((
"Adding set diffs for Col: %d.\n", col)) ;
1254 COLAMD_ASSERT(row >= 0 && row < n_row) ;
1255 row_mark = Row [row].shared2.mark ;
1257 if (ROW_IS_MARKED_DEAD (row_mark))
1261 COLAMD_ASSERT (row_mark > tag_mark) ;
1267 cur_score += row_mark - tag_mark ;
1269 cur_score = numext::mini(cur_score, n_col) ;
1273 Col [col].length = (IndexType) (new_cp - &A [Col [col].start]) ;
1277 if (Col [col].length == 0)
1279 COLAMD_DEBUG4 ((
"further mass elimination. Col: %d\n", col)) ;
1281 KILL_PRINCIPAL_COL (col) ;
1282 pivot_row_degree -= Col [col].shared1.thickness ;
1283 COLAMD_ASSERT (pivot_row_degree >= 0) ;
1285 Col [col].shared2.order = k ;
1287 k += Col [col].shared1.thickness ;
1293 COLAMD_DEBUG4 ((
"Preparing supercol detection for Col: %d.\n", col)) ;
1296 Col [col].shared2.score = cur_score ;
1301 COLAMD_DEBUG4 ((
" Hash = %d, n_col = %d.\n", hash, n_col)) ;
1302 COLAMD_ASSERT (hash <= n_col) ;
1304 head_column = head [hash] ;
1305 if (head_column > COLAMD_EMPTY)
1309 first_col = Col [head_column].shared3.headhash ;
1310 Col [head_column].shared3.headhash = col ;
1315 first_col = - (head_column + 2) ;
1316 head [hash] = - (col + 2) ;
1318 Col [col].shared4.hash_next = first_col ;
1321 Col [col].shared3.hash = (IndexType) hash ;
1322 COLAMD_ASSERT (COL_IS_ALIVE (col)) ;
1330 COLAMD_DEBUG3 ((
"** Supercolumn detection phase. **\n")) ;
1332 Eigen::internal::detect_super_cols (Col, A, head, pivot_row_start, pivot_row_length) ;
1336 KILL_PRINCIPAL_COL (pivot_col) ;
1340 tag_mark += (max_deg + 1) ;
1341 if (tag_mark >= max_mark)
1343 COLAMD_DEBUG2 ((
"clearing tag_mark\n")) ;
1344 tag_mark = Eigen::internal::clear_mark (n_row, Row) ;
1349 COLAMD_DEBUG3 ((
"** Finalize scores phase. **\n")) ;
1352 rp = &A [pivot_row_start] ;
1355 rp_end = rp + pivot_row_length ;
1360 if (COL_IS_DEAD (col))
1366 A [Col [col].start + (Col [col].length++)] = pivot_row ;
1371 cur_score = Col [col].shared2.score + pivot_row_degree ;
1376 max_score = n_col - k - Col [col].shared1.thickness ;
1379 cur_score -= Col [col].shared1.thickness ;
1382 cur_score = numext::mini(cur_score, max_score) ;
1383 COLAMD_ASSERT (cur_score >= 0) ;
1386 Col [col].shared2.score = cur_score ;
1390 COLAMD_ASSERT (min_score >= 0) ;
1391 COLAMD_ASSERT (min_score <= n_col) ;
1392 COLAMD_ASSERT (cur_score >= 0) ;
1393 COLAMD_ASSERT (cur_score <= n_col) ;
1394 COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ;
1395 next_col = head [cur_score] ;
1396 Col [col].shared4.degree_next = next_col ;
1397 Col [col].shared3.prev = COLAMD_EMPTY ;
1398 if (next_col != COLAMD_EMPTY)
1400 Col [next_col].shared3.prev = col ;
1402 head [cur_score] = col ;
1405 min_score = numext::mini(min_score, cur_score) ;
1411 if (pivot_row_degree > 0)
1415 Row [pivot_row].start = pivot_row_start ;
1416 Row [pivot_row].length = (IndexType) (new_rp - &A[pivot_row_start]) ;
1417 Row [pivot_row].shared1.degree = pivot_row_degree ;
1418 Row [pivot_row].shared2.mark = 0 ;
1445 template <
typename IndexType>
1446 static inline void order_children
1451 colamd_col<IndexType> Col [],
1464 for (i = 0 ; i < n_col ; i++)
1467 COLAMD_ASSERT (COL_IS_DEAD (i)) ;
1468 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY)
1474 parent = Col [parent].shared1.parent ;
1475 }
while (!COL_IS_DEAD_PRINCIPAL (parent)) ;
1481 order = Col [parent].shared2.order ;
1485 COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ;
1488 Col [c].shared2.order = order++ ;
1490 Col [c].shared1.parent = parent ;
1493 c = Col [c].shared1.parent ;
1498 }
while (Col [c].shared2.order == COLAMD_EMPTY) ;
1501 Col [parent].shared2.order = order ;
1507 for (c = 0 ; c < n_col ; c++)
1509 p [Col [c].shared2.order] = c ;
1546 template <
typename IndexType>
1547 static void detect_super_cols
1551 colamd_col<IndexType> Col [],
1554 IndexType row_start,
1555 IndexType row_length
1571 IndexType head_column ;
1572 IndexType first_col ;
1576 rp = &A [row_start] ;
1577 rp_end = rp + row_length ;
1581 if (COL_IS_DEAD (col))
1587 hash = Col [col].shared3.hash ;
1588 COLAMD_ASSERT (hash <= n_col) ;
1592 head_column = head [hash] ;
1593 if (head_column > COLAMD_EMPTY)
1595 first_col = Col [head_column].shared3.headhash ;
1599 first_col = - (head_column + 2) ;
1604 for (super_c = first_col ; super_c != COLAMD_EMPTY ;
1605 super_c = Col [super_c].shared4.hash_next)
1607 COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ;
1608 COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ;
1609 length = Col [super_c].length ;
1616 for (c = Col [super_c].shared4.hash_next ;
1617 c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next)
1619 COLAMD_ASSERT (c != super_c) ;
1620 COLAMD_ASSERT (COL_IS_ALIVE (c)) ;
1621 COLAMD_ASSERT (Col [c].shared3.hash == hash) ;
1624 if (Col [c].length != length ||
1625 Col [c].shared2.score != Col [super_c].shared2.score)
1632 cp1 = &A [Col [super_c].start] ;
1633 cp2 = &A [Col [c].start] ;
1635 for (i = 0 ; i < length ; i++)
1638 COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ;
1639 COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ;
1642 if (*cp1++ != *cp2++)
1657 COLAMD_ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
1659 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
1660 Col [c].shared1.parent = super_c ;
1661 KILL_NON_PRINCIPAL_COL (c) ;
1663 Col [c].shared2.order = COLAMD_EMPTY ;
1665 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
1671 if (head_column > COLAMD_EMPTY)
1674 Col [head_column].shared3.headhash = COLAMD_EMPTY ;
1679 head [hash] = COLAMD_EMPTY ;
1697 template <
typename IndexType>
1698 static IndexType garbage_collection
1704 Colamd_Row<IndexType> Row [],
1705 colamd_col<IndexType> Col [],
1722 for (c = 0 ; c < n_col ; c++)
1724 if (COL_IS_ALIVE (c))
1726 psrc = &A [Col [c].start] ;
1729 COLAMD_ASSERT (pdest <= psrc) ;
1730 Col [c].start = (IndexType) (pdest - &A [0]) ;
1731 length = Col [c].length ;
1732 for (j = 0 ; j < length ; j++)
1735 if (ROW_IS_ALIVE (r))
1740 Col [c].length = (IndexType) (pdest - &A [Col [c].start]) ;
1746 for (r = 0 ; r < n_row ; r++)
1748 if (ROW_IS_ALIVE (r))
1750 if (Row [r].length == 0)
1753 COLAMD_DEBUG3 ((
"Defrag row kill\n")) ;
1759 psrc = &A [Row [r].start] ;
1760 Row [r].shared2.first_column = *psrc ;
1761 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
1763 *psrc = ONES_COMPLEMENT (r) ;
1772 while (psrc < pfree)
1779 r = ONES_COMPLEMENT (*psrc) ;
1780 COLAMD_ASSERT (r >= 0 && r < n_row) ;
1782 *psrc = Row [r].shared2.first_column ;
1783 COLAMD_ASSERT (ROW_IS_ALIVE (r)) ;
1786 COLAMD_ASSERT (pdest <= psrc) ;
1787 Row [r].start = (IndexType) (pdest - &A [0]) ;
1788 length = Row [r].length ;
1789 for (j = 0 ; j < length ; j++)
1792 if (COL_IS_ALIVE (c))
1797 Row [r].length = (IndexType) (pdest - &A [Row [r].start]) ;
1802 COLAMD_ASSERT (debug_rows == 0) ;
1806 return ((IndexType) (pdest - &A [0])) ;
1818 template <
typename IndexType>
1819 static inline IndexType clear_mark
1824 Colamd_Row<IndexType> Row []
1831 for (r = 0 ; r < n_row ; r++)
1833 if (ROW_IS_ALIVE (r))
1835 Row [r].shared2.mark = 0 ;
Definition: Eigen_Colamd.h:50