47 #ifndef EIGEN_COLAMD_H
48 #define EIGEN_COLAMD_H
110 template <
typename IndexType>
135 template <
typename IndexType>
179 template <
typename IndexType>
221 template <
typename IndexType>
225 template <
typename IndexType>
230 template <
typename IndexType>
231 static IndexType
init_rows_cols (IndexType n_row, IndexType n_col, RowStructure<IndexType>
Row [], ColStructure<IndexType>
col [], IndexType
A [], IndexType
p [], IndexType
stats[
NStats] );
233 template <
typename IndexType>
234 static void init_scoring (IndexType n_row, IndexType n_col, RowStructure<IndexType>
Row [], ColStructure<IndexType> Col [], IndexType
A [], IndexType
head [],
double knobs[
NKnobs], IndexType *p_n_row2, IndexType *p_n_col2, IndexType *p_max_deg);
236 template <
typename IndexType>
237 static IndexType
find_ordering (IndexType n_row, IndexType n_col, IndexType Alen, RowStructure<IndexType>
Row [], ColStructure<IndexType> Col [], IndexType
A [], IndexType
head [], IndexType n_col2, IndexType max_deg, IndexType pfree);
239 template <
typename IndexType>
240 static void order_children (IndexType n_col, ColStructure<IndexType> Col [], IndexType
p []);
242 template <
typename IndexType>
243 static void detect_super_cols (ColStructure<IndexType> Col [], IndexType
A [], IndexType
head [], IndexType row_start, IndexType row_length ) ;
245 template <
typename IndexType>
246 static IndexType
garbage_collection (IndexType n_row, IndexType n_col, RowStructure<IndexType>
Row [], ColStructure<IndexType> Col [], IndexType
A [], IndexType *pfree) ;
248 template <
typename IndexType>
249 static inline IndexType
clear_mark (IndexType n_row, RowStructure<IndexType>
Row [] ) ;
253 #define COLAMD_DEBUG0(params) ;
254 #define COLAMD_DEBUG1(params) ;
255 #define COLAMD_DEBUG2(params) ;
256 #define COLAMD_DEBUG3(params) ;
257 #define COLAMD_DEBUG4(params) ;
259 #define COLAMD_ASSERT(expression) ((void) 0)
276 template <
typename IndexType>
277 inline IndexType
recommended ( IndexType nnz, IndexType n_row, IndexType n_col)
279 if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0)
282 return (2 * (nnz) +
colamd_c (n_col) +
colamd_r (n_row) + (n_col) + ((nnz) / 5));
341 template <
typename IndexType>
357 double default_knobs [
NKnobs] ;
410 COLAMD_DEBUG0 ((
"colamd: number of entries negative %d\n", nnz)) ;
427 knobs = default_knobs ;
434 need = 2*nnz + n_col + Col_size + Row_size ;
442 COLAMD_DEBUG0 ((
"colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
446 Alen -= Col_size + Row_size ;
462 &n_row2, &n_col2, &max_deg) ;
467 n_col2, max_deg, 2*nnz) ;
500 template <
typename IndexType>
531 if ((Col [
col].length) < 0)
541 Col [
col].shared1.thickness = 1 ;
542 Col [
col].shared2.score = 0 ;
544 Col [
col].shared4.degree_next =
Empty ;
556 Row [
row].shared2.mark = -1 ;
564 cp_end = &
A [
p [
col+1]] ;
571 if (row < 0 || row >= n_row)
614 Row [0].start =
p [n_col] ;
615 Row [0].shared1.p =
Row [0].start ;
616 Row [0].shared2.mark = -1 ;
621 Row [
row].shared2.mark = -1 ;
632 cp_end = &
A [
p [
col+1]] ;
650 cp_end = &
A [
p [
col+1]] ;
653 A [(
Row [*cp++].shared1.p)++] =
col ;
662 Row [
row].shared2.mark = 0 ;
670 COLAMD_DEBUG0 ((
"colamd: reconstructing column form, matrix jumbled\n")) ;
680 p [0] = Col [0].start ;
685 Col [
col].start = Col [
col-1].start + Col [
col-1].length ;
694 rp_end = rp +
Row [
row].length ;
697 A [(
p [*rp++])++] =
row ;
716 template <
typename IndexType>
741 IndexType col_length ;
745 IndexType dense_row_count ;
746 IndexType dense_col_count ;
747 IndexType min_score ;
756 COLAMD_DEBUG1 ((
"colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
765 for (
c = n_col-1 ;
c >= 0 ;
c--)
775 COLAMD_DEBUG1 ((
"colamd: null columns killed: %d\n", n_col - n_col2)) ;
780 for (
c = n_col-1 ;
c >= 0 ;
c--)
783 if (Col[
c].is_dead())
788 if (deg > dense_col_count)
794 cp_end = cp + Col [
c].length ;
797 Row [*cp++].shared1.degree-- ;
799 Col[
c].kill_principal() ;
802 COLAMD_DEBUG1 ((
"colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
806 for (r = 0 ; r < n_row ; r++)
808 deg =
Row [r].shared1.degree ;
810 if (deg > dense_row_count || deg == 0)
822 COLAMD_DEBUG1 ((
"colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
832 for (
c = n_col-1 ;
c >= 0 ;
c--)
835 if (Col[
c].is_dead())
842 cp_end = cp + Col [
c].length ;
855 score +=
Row [
row].shared1.degree - 1 ;
860 col_length = (IndexType) (new_cp - &
A [Col [
c].start]) ;
866 Col [
c].shared2.order = --n_col2 ;
867 Col[
c].kill_principal() ;
874 Col [
c].length = col_length ;
875 Col [
c].shared2.score = score ;
878 COLAMD_DEBUG1 ((
"colamd: Dense, null, and newly-null columns killed: %d\n",
890 for (
c = 0 ;
c <= n_col ;
c++)
897 for (
c = n_col-1 ;
c >= 0 ;
c--)
900 if (Col[
c].is_alive())
903 c, Col [
c].shared2.
score, min_score, n_col)) ;
916 next_col =
head [score] ;
922 if (next_col !=
Empty)
940 *p_max_deg = max_deg ;
953 template <
typename IndexType>
973 IndexType pivot_col ;
976 IndexType pivot_row ;
979 IndexType pivot_row_start ;
980 IndexType pivot_row_degree ;
981 IndexType pivot_row_length ;
982 IndexType pivot_col_score ;
983 IndexType needed_memory ;
988 IndexType max_score ;
989 IndexType cur_score ;
991 IndexType head_column ;
992 IndexType first_col ;
995 IndexType set_difference ;
996 IndexType min_score ;
997 IndexType col_thickness ;
999 IndexType pivot_col_thickness ;
1000 IndexType prev_col ;
1001 IndexType next_col ;
1002 IndexType ngarbage ;
1007 max_mark = INT_MAX - n_col ;
1015 for (k = 0 ; k < n_col2 ; )
1026 while (min_score < n_col &&
head [min_score] ==
Empty)
1030 pivot_col =
head [min_score] ;
1032 next_col = Col [pivot_col].shared4.degree_next ;
1033 head [min_score] = next_col ;
1034 if (next_col !=
Empty)
1036 Col [next_col].shared3.prev =
Empty ;
1043 pivot_col_score = Col [pivot_col].shared2.score ;
1046 Col [pivot_col].shared2.order = k ;
1049 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
1050 k += pivot_col_thickness ;
1055 needed_memory =
numext::mini(pivot_col_score, n_col - k) ;
1056 if (pfree + needed_memory >= Alen)
1070 pivot_row_start = pfree ;
1073 pivot_row_degree = 0 ;
1077 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
1080 cp = &
A [Col [pivot_col].start] ;
1081 cp_end = cp + Col [pivot_col].length ;
1093 rp_end = rp +
Row [
row].length ;
1099 col_thickness = Col [
col].shared1.thickness ;
1100 if (col_thickness > 0 && Col[
col].is_alive())
1103 Col [
col].shared1.thickness = -col_thickness ;
1107 pivot_row_degree += col_thickness ;
1113 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
1120 cp = &
A [Col [pivot_col].start] ;
1121 cp_end = cp + Col [pivot_col].length ;
1132 pivot_row_length = pfree - pivot_row_start ;
1133 if (pivot_row_length > 0)
1136 pivot_row =
A [Col [pivot_col].start] ;
1145 COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
1168 COLAMD_DEBUG3 ((
"** Computing set differences phase. **\n")) ;
1174 rp = &
A [pivot_row_start] ;
1175 rp_end = rp + pivot_row_length ;
1183 col_thickness = -Col [
col].shared1.thickness ;
1185 Col [
col].shared1.thickness = col_thickness ;
1189 cur_score = Col [
col].shared2.score ;
1190 prev_col = Col [
col].shared3.prev ;
1191 next_col = Col [
col].shared4.degree_next ;
1195 if (prev_col ==
Empty)
1197 head [cur_score] = next_col ;
1201 Col [prev_col].shared4.degree_next = next_col ;
1203 if (next_col !=
Empty)
1205 Col [next_col].shared3.prev = prev_col ;
1210 cp = &
A [Col [
col].start] ;
1211 cp_end = cp + Col [
col].length ;
1221 row_mark =
Row [
row].shared2.mark ;
1223 set_difference = row_mark - tag_mark ;
1225 if (set_difference < 0)
1228 set_difference =
Row [
row].shared1.degree ;
1231 set_difference -= col_thickness ;
1234 if (set_difference == 0)
1242 Row [
row].shared2.mark = set_difference + tag_mark ;
1253 rp = &
A [pivot_row_start] ;
1254 rp_end = rp + pivot_row_length ;
1262 cp = &
A [Col [
col].start] ;
1265 cp_end = cp + Col [
col].length ;
1279 row_mark =
Row [
row].shared2.mark ;
1286 cur_score += row_mark - tag_mark ;
1292 Col [
col].length = (IndexType) (new_cp - &
A [Col [
col].start]) ;
1296 if (Col [
col].length == 0)
1300 Col[
col].kill_principal() ;
1301 pivot_row_degree -= Col [
col].shared1.thickness ;
1304 Col [
col].shared2.order = k ;
1306 k += Col [
col].shared1.thickness ;
1315 Col [
col].shared2.score = cur_score ;
1324 if (head_column >
Empty)
1328 first_col = Col [head_column].shared3.headhash ;
1329 Col [head_column].shared3.headhash =
col ;
1334 first_col = - (head_column + 2) ;
1337 Col [
col].shared4.hash_next = first_col ;
1340 Col [
col].shared3.hash = (IndexType)
hash ;
1355 Col[pivot_col].kill_principal() ;
1359 tag_mark += (max_deg + 1) ;
1360 if (tag_mark >= max_mark)
1371 rp = &
A [pivot_row_start] ;
1374 rp_end = rp + pivot_row_length ;
1379 if (Col[
col].is_dead())
1385 A [Col [
col].start + (Col [
col].length++)] = pivot_row ;
1390 cur_score = Col [
col].shared2.score + pivot_row_degree ;
1395 max_score = n_col - k - Col [
col].shared1.thickness ;
1398 cur_score -= Col [
col].shared1.thickness ;
1405 Col [
col].shared2.score = cur_score ;
1414 next_col =
head [cur_score] ;
1415 Col [
col].shared4.degree_next = next_col ;
1417 if (next_col !=
Empty)
1419 Col [next_col].shared3.prev =
col ;
1430 if (pivot_row_degree > 0)
1434 Row [pivot_row].start = pivot_row_start ;
1435 Row [pivot_row].length = (IndexType) (new_rp - &
A[pivot_row_start]) ;
1436 Row [pivot_row].shared1.degree = pivot_row_degree ;
1437 Row [pivot_row].shared2.mark = 0 ;
1464 template <
typename IndexType>
1483 for (
i = 0 ;
i < n_col ;
i++)
1487 if (!Col[
i].is_dead_principal() && Col [
i].shared2.
order ==
Empty)
1494 }
while (!Col[parent].is_dead_principal()) ;
1526 for (
c = 0 ;
c < n_col ;
c++)
1565 template <
typename IndexType>
1573 IndexType row_start,
1574 IndexType row_length
1590 IndexType head_column ;
1591 IndexType first_col ;
1595 rp = &
A [row_start] ;
1596 rp_end = rp + row_length ;
1600 if (Col[
col].is_dead())
1612 if (head_column >
Empty)
1618 first_col = - (head_column + 2) ;
1623 for (super_c = first_col ; super_c !=
Empty ;
1628 length = Col [super_c].
length ;
1643 if (Col [
c].length != length ||
1651 cp1 = &
A [Col [super_c].
start] ;
1652 cp2 = &
A [Col [
c].start] ;
1654 for (
i = 0 ;
i < length ;
i++)
1661 if (*cp1++ != *cp2++)
1676 COLAMD_ASSERT (Col [
c].shared2.score == Col [super_c].shared2.score) ;
1678 Col [super_c].shared1.thickness += Col [
c].shared1.thickness ;
1679 Col [
c].shared1.parent = super_c ;
1680 Col[
c].kill_non_principal() ;
1682 Col [
c].shared2.order =
Empty ;
1684 Col [prev_c].shared4.hash_next = Col [
c].shared4.hash_next ;
1690 if (head_column >
Empty)
1716 template <
typename IndexType>
1741 for (
c = 0 ;
c < n_col ;
c++)
1743 if (Col[
c].is_alive())
1745 psrc = &
A [Col [
c].start] ;
1749 Col [
c].start = (IndexType) (pdest - &
A [0]) ;
1750 length = Col [
c].length ;
1751 for (
j = 0 ;
j < length ;
j++)
1754 if (
Row[r].is_alive())
1759 Col [
c].length = (IndexType) (pdest - &
A [Col [
c].start]) ;
1765 for (r = 0 ; r < n_row ; r++)
1767 if (
Row[r].is_alive())
1769 if (
Row [r].length == 0)
1778 psrc = &
A [
Row [r].start] ;
1779 Row [r].shared2.first_column = *psrc ;
1791 while (psrc < pfree)
1801 *psrc =
Row [r].shared2.first_column ;
1806 Row [r].start = (IndexType) (pdest - &
A [0]) ;
1807 length =
Row [r].length ;
1808 for (
j = 0 ;
j < length ;
j++)
1811 if (Col[
c].is_alive())
1816 Row [r].length = (IndexType) (pdest - &
A [
Row [r].start]) ;
1825 return ((IndexType) (pdest - &
A [0])) ;
1837 template <
typename IndexType>
1850 for (r = 0 ; r < n_row ; r++)
1852 if (
Row[r].is_alive())
1854 Row [r].shared2.mark = 0 ;