733 #ifdef MATLAB_MEX_FILE
747 #define PRIVATE static
749 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
750 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
752 #define ONES_COMPLEMENT(r) (-(r)-1)
775 #define DEAD_PRINCIPAL (-1)
776 #define DEAD_NON_PRINCIPAL (-2)
779 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
780 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
781 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
782 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
783 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
784 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
785 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
786 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
787 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
793 #ifdef MATLAB_MEX_FILE
796 #define PRINTF mexPrintf
799 #define INDEX(i) ((i)+1)
806 #define PRINTF printf
911 #define DEBUG0(params) { (void) PRINTF params ; }
912 #define DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; }
913 #define DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; }
914 #define DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; }
915 #define DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; }
917 #ifdef MATLAB_MEX_FILE
918 #define ASSERT(expression) (mxAssert ((expression), ""))
920 #define ASSERT(expression) (assert (expression))
971 #define DEBUG0(params) ;
972 #define DEBUG1(params) ;
973 #define DEBUG2(params) ;
974 #define DEBUG3(params) ;
975 #define DEBUG4(params) ;
977 #define ASSERT(expression) ((void) 0)
1075 void * (*allocate) (size_t, size_t),
1078 void (*release) (
void *)
1104 colamd_get_debug (
"symamd") ;
1111 DEBUG0 ((
"symamd: stats not present\n")) ;
1125 DEBUG0 ((
"symamd: A not present\n")) ;
1132 DEBUG0 ((
"symamd: p not present\n")) ;
1140 DEBUG0 ((
"symamd: n negative %d\n", n)) ;
1149 DEBUG0 ((
"symamd: number of entries negative %d\n", nnz)) ;
1157 DEBUG0 ((
"symamd: p[0] not zero %d\n", p [0])) ;
1166 knobs = default_knobs ;
1171 count = (
integer *) ((*allocate) (n+1,
sizeof (
int))) ;
1175 DEBUG0 ((
"symamd: allocate count (size %d) failed\n", n+1)) ;
1179 mark = (
integer *) ((*allocate) (n+1,
sizeof (
int))) ;
1183 (*release) ((
void *) count) ;
1184 DEBUG0 ((
"symamd: allocate mark (size %d) failed\n", n+1)) ;
1192 for (i = 0 ; i < n ; i++)
1197 for (j = 0 ; j < n ; j++)
1201 length = p [j+1] - p [j] ;
1208 (*release) ((
void *) count) ;
1209 (*release) ((
void *) mark) ;
1210 DEBUG0 ((
"symamd: col %d negative length %d\n", j, length)) ;
1214 for (pp = p [j] ; pp < p [j+1] ; pp++)
1217 if (i < 0 || i >= n)
1224 (*release) ((
void *) count) ;
1225 (*release) ((
void *) mark) ;
1226 DEBUG0 ((
"symamd: row %d col %d out of bounds\n", i, j)) ;
1230 if (i <= last_row || mark [i] == j)
1238 DEBUG1 ((
"symamd: row %d col %d unsorted/duplicate\n", i, j)) ;
1241 if (i > j && mark [i] != j)
1258 (*release) ((
void *) mark) ;
1265 for (j = 1 ; j <= n ; j++)
1267 perm [j] = perm [j-1] + count [j-1] ;
1269 for (j = 0 ; j < n ; j++)
1271 count [j] = perm [j] ;
1279 M = (
integer *) ((*allocate) (Mlen,
sizeof (
int))) ;
1280 DEBUG0 ((
"symamd: M is %d-by-%d with %d entries, Mlen = %d\n",
1281 n_row, n, mnz, Mlen)) ;
1286 (*release) ((
void *) count) ;
1287 (*release) ((
void *) mark) ;
1288 DEBUG0 ((
"symamd: allocate M (size %d) failed\n", Mlen)) ;
1297 for (j = 0 ; j < n ; j++)
1299 ASSERT (p [j+1] - p [j] >= 0) ;
1300 for (pp = p [j] ; pp < p [j+1] ; pp++)
1303 ASSERT (i >= 0 && i < n) ;
1307 M [count [i]++] = k ;
1308 M [count [j]++] = k ;
1317 DEBUG0 ((
"symamd: Duplicates in A.\n")) ;
1318 for (i = 0 ; i < n ; i++)
1322 for (j = 0 ; j < n ; j++)
1324 ASSERT (p [j+1] - p [j] >= 0) ;
1325 for (pp = p [j] ; pp < p [j+1] ; pp++)
1328 ASSERT (i >= 0 && i < n) ;
1329 if (i > j && mark [i] != j)
1332 M [count [i]++] = k ;
1333 M [count [j]++] = k ;
1339 (*release) ((
void *) mark) ;
1343 (*release) ((
void *) count) ;
1350 cknobs [i] = knobs [i] ;
1356 if (n_row != 0 && n < n_row)
1372 if (!
mbdyn_colamd (n_row, n, Mlen, M, perm, cknobs, cstats))
1376 (*release) ((
void *) M) ;
1377 DEBUG0 ((
"symamd: internal error!\n")) ;
1392 (*release) ((
void *) M) ;
1393 DEBUG0 ((
"symamd: done.\n")) ;
1439 colamd_get_debug (
"colamd") ;
1446 DEBUG0 ((
"colamd: stats not present\n")) ;
1460 DEBUG0 ((
"colamd: A not present\n")) ;
1467 DEBUG0 ((
"colamd: p not present\n")) ;
1475 DEBUG0 ((
"colamd: nrow negative %d\n", n_row)) ;
1483 DEBUG0 ((
"colamd: ncol negative %d\n", n_col)) ;
1492 DEBUG0 ((
"colamd: number of entries negative %d\n", nnz)) ;
1500 DEBUG0 ((
"colamd: p[0] not zero %d\n", p [0])) ;
1509 knobs = default_knobs ;
1516 need = 2*nnz + n_col + Col_size + Row_size ;
1524 DEBUG0 ((
"colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen));
1528 Alen -= Col_size + Row_size ;
1537 DEBUG0 ((
"colamd: Matrix invalid\n")) ;
1544 &n_row2, &n_col2, &max_deg) ;
1548 ngarbage =
find_ordering (n_row, n_col, Alen, Row, Col, A, p,
1549 n_col2, max_deg, 2*nnz) ;
1560 DEBUG0 ((
"colamd: done.\n")) ;
1637 for (col = 0 ; col < n_col ; col++)
1639 Col [col].start = p [col] ;
1640 Col [col].length = p [col+1] - p [col] ;
1642 if (Col [col].length < 0)
1648 DEBUG0 ((
"colamd: col %d length %d < 0\n", col, Col [col].length)) ;
1652 Col [col].shared1.thickness = 1 ;
1653 Col [col].shared2.score = 0 ;
1654 Col [col].shared3.prev =
EMPTY ;
1655 Col [col].shared4.degree_next =
EMPTY ;
1664 for (row = 0 ; row < n_row ; row++)
1666 Row [row].length = 0 ;
1667 Row [row].shared2.mark = -1 ;
1670 for (col = 0 ; col < n_col ; col++)
1675 cp_end = &A [p [col+1]] ;
1682 if (row < 0 || row >= n_row)
1688 DEBUG0 ((
"colamd: row %d col %d out of bounds\n", row, col)) ;
1692 if (row <= last_row || Row [row].shared2.mark == col)
1700 DEBUG1 ((
"colamd: row %d col %d unsorted/duplicate\n",row,col));
1703 if (Row [row].shared2.mark != col)
1705 Row [row].length++ ;
1711 Col [col].length-- ;
1715 Row [row].shared2.mark = col ;
1725 Row [0].start = p [n_col] ;
1726 Row [0].shared1.p = Row [0].start ;
1727 Row [0].shared2.mark = -1 ;
1728 for (row = 1 ; row < n_row ; row++)
1730 Row [row].start = Row [row-1].start + Row [row-1].length ;
1731 Row [row].shared1.p = Row [row].start ;
1732 Row [row].shared2.mark = -1 ;
1740 for (col = 0 ; col < n_col ; col++)
1743 cp_end = &A [p [col+1]] ;
1747 if (Row [row].shared2.mark != col)
1749 A [(Row [row].shared1.p)++] = col ;
1750 Row [row].shared2.mark = col ;
1758 for (col = 0 ; col < n_col ; col++)
1761 cp_end = &A [p [col+1]] ;
1764 A [(Row [*cp++].shared1.p)++] = col ;
1771 for (row = 0 ; row < n_row ; row++)
1773 Row [row].shared2.mark = 0 ;
1774 Row [row].shared1.degree = Row [row].length ;
1781 DEBUG0 ((
"colamd: reconstructing column form, matrix jumbled\n")) ;
1785 for (col = 0 ; col < n_col ; col++)
1787 p [col] = Col [col].length ;
1789 for (row = 0 ; row < n_row ; row++)
1791 rp = &A [Row [row].start] ;
1792 rp_end = rp + Row [row].length ;
1798 for (col = 0 ; col < n_col ; col++)
1812 p [0] = Col [0].start ;
1813 for (col = 1 ; col < n_col ; col++)
1817 Col [col].start = Col [col-1].start + Col [col-1].length ;
1818 p [col] = Col [col].start ;
1823 for (row = 0 ; row < n_row ; row++)
1825 rp = &A [Row [row].start] ;
1826 rp_end = rp + Row [row].length ;
1829 A [(p [*rp++])++] = row ;
1859 double knobs [COLAMD_KNOBS],
1891 DEBUG1 ((
"colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ;
1900 for (c = n_col-1 ; c >= 0 ; c--)
1910 DEBUG1 ((
"colamd: null columns killed: %d\n", n_col - n_col2)) ;
1915 for (c = n_col-1 ; c >= 0 ; c--)
1923 if (deg > dense_col_count)
1928 cp = &A [Col [
c].
start] ;
1929 cp_end = cp + Col [
c].length ;
1937 DEBUG1 ((
"colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ;
1941 for (r = 0 ; r < n_row ; r++)
1944 ASSERT (deg >= 0 && deg <= n_col) ;
1945 if (deg > dense_row_count || deg == 0)
1954 max_deg =
MAX (max_deg, deg) ;
1957 DEBUG1 ((
"colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ;
1967 for (c = n_col-1 ; c >= 0 ; c--)
1975 cp = &A [Col [
c].
start] ;
1977 cp_end = cp + Col [
c].length ;
1992 score =
MIN (score, n_col) ;
1995 col_length = (
int) (new_cp - &A [Col [c].start]) ;
1996 if (col_length == 0)
2000 DEBUG2 ((
"Newly null killed: %d\n", c)) ;
2001 Col [
c].shared2.order = --n_col2 ;
2008 ASSERT (score <= n_col) ;
2009 Col [
c].length = col_length ;
2010 Col [
c].shared2.score = score ;
2013 DEBUG1 ((
"colamd: Dense, null, and newly-null columns killed: %d\n",
2022 debug_structures (n_row, n_col, Row, Col, A, n_col2) ;
2032 for (c = 0 ; c <= n_col ; c++)
2039 for (c = n_col-1 ; c >= 0 ; c--)
2044 DEBUG4 ((
"place %d score %d minscore %d ncol %d\n",
2045 c, Col [c].shared2.
score, min_score, n_col)) ;
2051 ASSERT (min_score >= 0) ;
2052 ASSERT (min_score <= n_col) ;
2054 ASSERT (score <= n_col) ;
2058 next_col = head [score] ;
2064 if (next_col !=
EMPTY)
2071 min_score =
MIN (min_score, score) ;
2081 DEBUG1 ((
"colamd: Live cols %d out of %d, non-princ: %d\n",
2082 debug_count, n_col, n_col-debug_count)) ;
2083 ASSERT (debug_count == n_col2) ;
2084 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ;
2089 *p_n_col2 = n_col2 ;
2090 *p_n_row2 = n_row2 ;
2091 *p_max_deg = max_deg ;
2162 max_mark = INT_MAX - n_col ;
2166 DEBUG1 ((
"colamd: Ordering, n_col2=%d\n", n_col2)) ;
2170 for (k = 0 ; k < n_col2 ; )
2174 if (debug_step % 100 == 0)
2176 DEBUG2 ((
"\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2180 DEBUG3 ((
"\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ;
2183 debug_deg_lists (n_row, n_col, Row, Col, head,
2184 min_score, n_col2-k, max_deg) ;
2185 debug_matrix (n_row, n_col, Row, Col, A) ;
2191 ASSERT (min_score >= 0) ;
2192 ASSERT (min_score <= n_col) ;
2196 for (debug_d = 0 ; debug_d < min_score ; debug_d++)
2203 while (head [min_score] ==
EMPTY && min_score < n_col)
2207 pivot_col = head [min_score] ;
2208 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
2209 next_col = Col [pivot_col].shared4.degree_next ;
2210 head [min_score] = next_col ;
2211 if (next_col !=
EMPTY)
2213 Col [next_col].shared3.prev =
EMPTY ;
2217 DEBUG3 ((
"Pivot col: %d\n", pivot_col)) ;
2220 pivot_col_score = Col [pivot_col].shared2.score ;
2223 Col [pivot_col].shared2.order = k ;
2226 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
2227 k += pivot_col_thickness ;
2228 ASSERT (pivot_col_thickness > 0) ;
2232 needed_memory =
MIN (pivot_col_score, n_col - k) ;
2233 if (pfree + needed_memory >= Alen)
2238 ASSERT (pfree + needed_memory < Alen) ;
2243 debug_matrix (n_row, n_col, Row, Col, A) ;
2250 pivot_row_start = pfree ;
2253 pivot_row_degree = 0 ;
2257 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
2260 cp = &A [Col [pivot_col].start] ;
2261 cp_end = cp + Col [pivot_col].length ;
2272 rp = &A [Row [row].start] ;
2273 rp_end = rp + Row [row].length ;
2279 col_thickness = Col [col].shared1.thickness ;
2283 Col [col].shared1.thickness = -col_thickness ;
2287 pivot_row_degree += col_thickness ;
2293 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
2294 max_deg =
MAX (max_deg, pivot_row_degree) ;
2298 debug_mark (n_row, Row, tag_mark, max_mark) ;
2304 cp = &A [Col [pivot_col].start] ;
2305 cp_end = cp + Col [pivot_col].length ;
2310 DEBUG3 ((
"Kill row in pivot col: %d\n", row)) ;
2316 pivot_row_length = pfree - pivot_row_start ;
2317 if (pivot_row_length > 0)
2320 pivot_row = A [Col [pivot_col].start] ;
2321 DEBUG3 ((
"Pivotal row is %d\n", pivot_row)) ;
2327 ASSERT (pivot_row_length == 0) ;
2329 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
2352 DEBUG3 ((
"** Computing set differences phase. **\n")) ;
2356 DEBUG3 ((
"Pivot row: ")) ;
2358 rp = &A [pivot_row_start] ;
2359 rp_end = rp + pivot_row_length ;
2364 DEBUG3 ((
"Col: %d\n", col)) ;
2367 col_thickness = -Col [col].shared1.thickness ;
2368 ASSERT (col_thickness > 0) ;
2369 Col [col].shared1.thickness = col_thickness ;
2373 cur_score = Col [col].shared2.score ;
2374 prev_col = Col [col].shared3.prev ;
2375 next_col = Col [col].shared4.degree_next ;
2376 ASSERT (cur_score >= 0) ;
2377 ASSERT (cur_score <= n_col) ;
2379 if (prev_col ==
EMPTY)
2381 head [cur_score] = next_col ;
2385 Col [prev_col].shared4.degree_next = next_col ;
2387 if (next_col !=
EMPTY)
2389 Col [next_col].shared3.prev = prev_col ;
2394 cp = &A [Col [col].start] ;
2395 cp_end = cp + Col [col].length ;
2400 row_mark = Row [row].shared2.mark ;
2406 ASSERT (row != pivot_row) ;
2407 set_difference = row_mark - tag_mark ;
2409 if (set_difference < 0)
2411 ASSERT (Row [row].shared1.degree <= max_deg) ;
2412 set_difference = Row [row].shared1.degree ;
2415 set_difference -= col_thickness ;
2416 ASSERT (set_difference >= 0) ;
2418 if (set_difference == 0)
2420 DEBUG3 ((
"aggressive absorption. Row: %d\n", row)) ;
2426 Row [row].shared2.mark = set_difference + tag_mark ;
2432 debug_deg_lists (n_row, n_col, Row, Col, head,
2433 min_score, n_col2-k-pivot_row_degree, max_deg) ;
2438 DEBUG3 ((
"** Adding set differences phase. **\n")) ;
2441 rp = &A [pivot_row_start] ;
2442 rp_end = rp + pivot_row_length ;
2450 cp = &A [Col [col].start] ;
2453 cp_end = cp + Col [col].length ;
2455 DEBUG4 ((
"Adding set diffs for Col: %d.\n", col)) ;
2461 ASSERT(row >= 0 && row < n_row) ;
2462 row_mark = Row [row].shared2.mark ;
2468 ASSERT (row_mark > tag_mark) ;
2474 cur_score += row_mark - tag_mark ;
2476 cur_score =
MIN (cur_score, n_col) ;
2480 Col [col].length = (
int) (new_cp - &A [Col [col].start]) ;
2484 if (Col [col].length == 0)
2486 DEBUG4 ((
"further mass elimination. Col: %d\n", col)) ;
2489 pivot_row_degree -= Col [col].shared1.thickness ;
2490 ASSERT (pivot_row_degree >= 0) ;
2492 Col [col].shared2.order = k ;
2494 k += Col [col].shared1.thickness ;
2500 DEBUG4 ((
"Preparing supercol detection for Col: %d.\n", col)) ;
2503 Col [col].shared2.score = cur_score ;
2508 DEBUG4 ((
" Hash = %d, n_col = %d.\n", hash, n_col)) ;
2511 head_column = head [hash] ;
2512 if (head_column >
EMPTY)
2516 first_col = Col [head_column].shared3.headhash ;
2517 Col [head_column].shared3.headhash = col ;
2522 first_col = - (head_column + 2) ;
2523 head [hash] = - (col + 2) ;
2525 Col [col].shared4.hash_next = first_col ;
2528 Col [col].shared3.hash = (
int) hash ;
2537 DEBUG3 ((
"** Supercolumn detection phase. **\n")) ;
2545 Col, A, head, pivot_row_start, pivot_row_length) ;
2553 tag_mark += (max_deg + 1) ;
2554 if (tag_mark >= max_mark)
2556 DEBUG2 ((
"clearing tag_mark\n")) ;
2562 debug_mark (n_row, Row, tag_mark, max_mark) ;
2567 DEBUG3 ((
"** Finalize scores phase. **\n")) ;
2570 rp = &A [pivot_row_start] ;
2573 rp_end = rp + pivot_row_length ;
2584 A [Col [col].start + (Col [col].length++)] = pivot_row ;
2589 cur_score = Col [col].shared2.score + pivot_row_degree ;
2594 max_score = n_col - k - Col [col].shared1.thickness ;
2597 cur_score -= Col [col].shared1.thickness ;
2600 cur_score =
MIN (cur_score, max_score) ;
2601 ASSERT (cur_score >= 0) ;
2604 Col [col].shared2.score = cur_score ;
2608 ASSERT (min_score >= 0) ;
2609 ASSERT (min_score <= n_col) ;
2610 ASSERT (cur_score >= 0) ;
2611 ASSERT (cur_score <= n_col) ;
2613 next_col = head [cur_score] ;
2614 Col [col].shared4.degree_next = next_col ;
2615 Col [col].shared3.prev =
EMPTY ;
2616 if (next_col !=
EMPTY)
2618 Col [next_col].shared3.prev = col ;
2620 head [cur_score] = col ;
2623 min_score =
MIN (min_score, cur_score) ;
2628 debug_deg_lists (n_row, n_col, Row, Col, head,
2629 min_score, n_col2-k, max_deg) ;
2634 if (pivot_row_degree > 0)
2638 Row [pivot_row].start = pivot_row_start ;
2639 Row [pivot_row].length = (
int) (new_rp - &A[pivot_row_start]) ;
2640 Row [pivot_row].shared1.degree = pivot_row_degree ;
2641 Row [pivot_row].shared2.mark = 0 ;
2687 for (i = 0 ; i < n_col ; i++)
2730 for (c = 0 ; c < n_col ; c++)
2805 rp = &A [row_start] ;
2806 rp_end = rp + row_length ;
2816 hash = Col [col].shared3.hash ;
2821 head_column = head [hash] ;
2822 if (head_column >
EMPTY)
2824 first_col = Col [head_column].shared3.headhash ;
2828 first_col = - (head_column + 2) ;
2833 for (super_c = first_col ; super_c !=
EMPTY ;
2834 super_c = Col [super_c].shared4.hash_next)
2837 ASSERT (Col [super_c].shared3.hash == hash) ;
2838 length = Col [super_c].length ;
2845 for (c = Col [super_c].shared4.hash_next ;
2846 c !=
EMPTY ; c = Col [c].shared4.hash_next)
2850 ASSERT (Col [c].shared3.hash == hash) ;
2853 if (Col [c].length != length ||
2854 Col [c].shared2.score != Col [super_c].shared2.score)
2861 cp1 = &A [Col [super_c].start] ;
2862 cp2 = &A [Col [
c].start] ;
2864 for (i = 0 ; i < length ; i++)
2871 if (*cp1++ != *cp2++)
2886 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
2888 Col [super_c].shared1.thickness += Col [
c].shared1.thickness ;
2889 Col [
c].shared1.parent = super_c ;
2892 Col [
c].shared2.order =
EMPTY ;
2894 Col [prev_c].shared4.hash_next = Col [
c].shared4.hash_next ;
2900 if (head_column >
EMPTY)
2903 Col [head_column].shared3.headhash =
EMPTY ;
2908 head [hash] =
EMPTY ;
2950 DEBUG2 ((
"Defrag..\n")) ;
2951 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
2958 for (c = 0 ; c < n_col ; c++)
2962 psrc = &A [Col [
c].start] ;
2966 Col [
c].start = (
int) (pdest - &A [0]) ;
2967 length = Col [
c].length ;
2968 for (j = 0 ; j < length ; j++)
2976 Col [
c].length = (
int) (pdest - &A [Col [c].start]) ;
2982 for (r = 0 ; r < n_row ; r++)
2986 if (Row [r].length == 0)
2989 DEBUG3 ((
"Defrag row kill\n")) ;
2995 psrc = &A [Row [r].start] ;
2996 Row [r].shared2.first_column = *psrc ;
3012 while (psrc < pfree)
3020 ASSERT (r >= 0 && r < n_row) ;
3022 *psrc = Row [r].shared2.first_column ;
3027 Row [r].start = (
int) (pdest - &A [0]) ;
3028 length = Row [r].length ;
3029 for (j = 0 ; j < length ; j++)
3037 Row [r].length = (
int) (pdest - &A [Row [r].start]) ;
3046 ASSERT (debug_rows == 0) ;
3050 return ((
int) (pdest - &A [0])) ;
3075 for (r = 0 ; r < n_row ; r++)
3079 Row [r].shared2.mark = 0 ;
3101 PRINTF (
"%s: No statistics available.\n", method) ;
3111 PRINTF (
"%s: OK. ", method) ;
3115 PRINTF (
"%s: ERROR. ", method) ;
3123 PRINTF (
"Matrix has unsorted or duplicate row indices.\n") ;
3125 PRINTF (
"%s: number of duplicate or out-of-order row indices: %d\n",
3128 PRINTF (
"%s: last seen duplicate or out-of-order row index: %d\n",
3129 method,
INDEX (i2)) ;
3131 PRINTF (
"%s: last seen in column: %d",
3132 method,
INDEX (i1)) ;
3140 PRINTF (
"%s: number of dense or empty rows ignored: %d\n",
3143 PRINTF (
"%s: number of dense or empty columns ignored: %d\n",
3146 PRINTF (
"%s: number of garbage collections performed: %d\n",
3152 PRINTF (
"Array A (row indices of matrix) not present.\n") ;
3157 PRINTF (
"Array p (column pointers for matrix) not present.\n") ;
3162 PRINTF (
"Invalid number of rows (%d).\n", i1) ;
3167 PRINTF (
"Invalid number of columns (%d).\n", i1) ;
3172 PRINTF (
"Invalid number of nonzero entries (%d).\n", i1) ;
3177 PRINTF (
"Invalid column pointer, p [0] = %d, must be zero.\n", i1) ;
3182 PRINTF (
"Array A too small.\n") ;
3183 PRINTF (
" Need Alen >= %d, but given only Alen = %d.\n",
3190 (
"Column %d has a negative number of nonzero entries (%d).\n",
3197 (
"Row index (row %d) out of bounds (%d to %d) in column %d.\n",
3203 PRINTF (
"Out of memory.\n") ;
3210 (
"Internal error! Please contact authors (davis@cise.ufl.edu).\n") ;
3265 for (c = 0 ; c < n_col ; c++)
3271 DEBUG4 ((
"initial live col %5d %5d %5d\n", c, len, score)) ;
3275 cp = &A [Col [
c].
start] ;
3285 i = Col [
c].shared2.order ;
3286 ASSERT (i >= n_col2 && i < n_col) ;
3290 for (r = 0 ; r < n_row ; r++)
3299 rp = &A [Row [r].
start] ;
3348 if (n_col > 10000 && colamd_debug <= 0)
3353 DEBUG4 ((
"Degree lists: %d\n", min_score)) ;
3354 for (deg = 0 ; deg <= n_col ; deg++)
3362 while (col !=
EMPTY)
3371 DEBUG4 ((
"should %d have %d\n", should, have)) ;
3372 ASSERT (should == have) ;
3376 if (n_row > 10000 && colamd_debug <= 0)
3380 for (row = 0 ; row < n_row ; row++)
3415 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
3416 if (n_row > 10000 && colamd_debug <= 0)
3420 for (r = 0 ; r < n_row ; r++)
3457 if (colamd_debug < 3)
3461 DEBUG3 ((
"DUMP MATRIX:\n")) ;
3462 for (r = 0 ; r < n_row ; r++)
3469 DEBUG3 ((
"start %d length %d degree %d\n",
3470 Row [r].start, Row [r].length, Row [r].shared1.
degree)) ;
3471 rp = &A [Row [r].
start] ;
3472 rp_end = rp + Row [r].length ;
3480 for (c = 0 ; c < n_col ; c++)
3487 DEBUG3 ((
"start %d length %d shared1 %d shared2 %d\n",
3488 Col [c].start, Col [c].length,
3490 cp = &A [Col [
c].
start] ;
3491 cp_end = cp + Col [
c].length ;
3510 colamd_debug = atoi (getenv (
"D")) ;
3513 DEBUG0 ((
"%s: debug version, D = %d (THIS WILL BE SLOW!)\n",
3514 method, colamd_debug)) ;
#define COLAMD_DEFRAG_COUNT
static integer clear_mark(integer n_row, mbdyn_Colamd_Row Row[])
static void order_children(integer n_col, mbdyn_Colamd_Col Col[], integer p[])
static integer init_rows_cols(integer n_row, integer n_col, mbdyn_Colamd_Row Row[], mbdyn_Colamd_Col Col[], integer A[], integer p[], integer stats[20])
static void init_scoring(integer n_row, integer n_col, mbdyn_Colamd_Row Row[], mbdyn_Colamd_Col Col[], integer A[], integer head[], double knobs[20], integer *p_n_row2, integer *p_n_col2, integer *p_max_deg)
#define COLAMD_ERROR_col_length_negative
union mbdyn_Colamd_Col::@1 shared2
#define COLAMD_RECOMMENDED(nnz, n_row, n_col)
void mbdyn_colamd_set_defaults(double knobs[20])
#define KILL_NON_PRINCIPAL_COL(c)
static integer find_ordering(integer n_row, integer n_col, integer Alen, mbdyn_Colamd_Row Row[], mbdyn_Colamd_Col Col[], integer A[], integer head[], integer n_col2, integer max_deg, integer pfree)
#define COLAMD_ERROR_p_not_present
#define COLAMD_ERROR_A_too_small
#define COLAMD_ERROR_A_not_present
void mbdyn_symamd_report(integer stats[20])
void mbdyn_colamd_report(integer stats[20])
#define COLAMD_ERROR_internal_error
integer mbdyn_colamd_recommended(integer nnz, integer n_row, integer n_col)
integer mbdyn_colamd(integer n_row, integer n_col, integer Alen, integer A[], integer p[], double knobs[20], integer stats[20])
union mbdyn_Colamd_Row::@4 shared1
union mbdyn_Colamd_Col::@0 shared1
#define COL_IS_DEAD_PRINCIPAL(c)
static void detect_super_cols(mbdyn_Colamd_Col Col[], integer A[], integer head[], integer row_start, integer row_length)
#define COLAMD_OK_BUT_JUMBLED
#define ASSERT(expression)
static std::stack< cleanup * > c
integer mbdyn_symamd(integer n, integer A[], integer p[], integer perm[], double knobs[20], integer stats[20], void *(*allocate)(size_t, size_t), void(*release)(void *))
union mbdyn_Colamd_Col::@3 shared4
#define KILL_PRINCIPAL_COL(c)
#define COLAMD_ERROR_ncol_negative
#define COLAMD_ERROR_row_index_out_of_bounds
#define ONES_COMPLEMENT(r)
#define COLAMD_ERROR_nrow_negative
union mbdyn_Colamd_Col::@2 shared3
#define COLAMD_ERROR_out_of_memory
static integer garbage_collection(integer n_row, integer n_col, mbdyn_Colamd_Row Row[], mbdyn_Colamd_Col Col[], integer A[], integer *pfree)
static void print_report(char *method, integer stats[20])
#define ROW_IS_MARKED_DEAD(row_mark)
#define COLAMD_ERROR_nnz_negative
#define COLAMD_ERROR_p0_nonzero