|
5164
|
1 /* ========================================================================== */ |
|
|
2 /* === colamd/symamd - a sparse matrix column ordering algorithm ============ */ |
|
|
3 /* ========================================================================== */ |
|
|
4 |
|
|
5 /* |
|
|
6 colamd: an approximate minimum degree column ordering algorithm, |
|
|
7 for LU factorization of symmetric or unsymmetric matrices, |
|
|
8 QR factorization, least squares, interior point methods for |
|
|
9 linear programming problems, and other related problems. |
|
|
10 |
|
|
11 symamd: an approximate minimum degree ordering algorithm for Cholesky |
|
|
12 factorization of symmetric matrices. |
|
|
13 |
|
|
14 Purpose: |
|
|
15 |
|
|
16 Colamd computes a permutation Q such that the Cholesky factorization of |
|
|
17 (AQ)'(AQ) has less fill-in and requires fewer floating point operations |
|
|
18 than A'A. This also provides a good ordering for sparse partial |
|
|
19 pivoting methods, P(AQ) = LU, where Q is computed prior to numerical |
|
|
20 factorization, and P is computed during numerical factorization via |
|
|
21 conventional partial pivoting with row interchanges. Colamd is the |
|
|
22 column ordering method used in SuperLU, part of the ScaLAPACK library. |
|
|
23 It is also available as built-in function in MATLAB Version 6, |
|
|
24 available from MathWorks, Inc. (http://www.mathworks.com). This |
|
|
25 routine can be used in place of colmmd in MATLAB. |
|
|
26 |
|
|
27 Symamd computes a permutation P of a symmetric matrix A such that the |
|
|
28 Cholesky factorization of PAP' has less fill-in and requires fewer |
|
|
29 floating point operations than A. Symamd constructs a matrix M such |
|
|
30 that M'M has the same nonzero pattern of A, and then orders the columns |
|
|
31 of M using colmmd. The column ordering of M is then returned as the |
|
|
32 row and column ordering P of A. |
|
|
33 |
|
|
34 Authors: |
|
|
35 |
|
|
36 The authors of the code itself are Stefan I. Larimore and Timothy A. |
|
|
37 Davis (davis@cise.ufl.edu), University of Florida. The algorithm was |
|
|
38 developed in collaboration with John Gilbert, Xerox PARC, and Esmond |
|
|
39 Ng, Oak Ridge National Laboratory. |
|
|
40 |
|
|
41 Date: |
|
|
42 |
|
|
43 September 8, 2003. Version 2.3. |
|
|
44 |
|
|
45 Acknowledgements: |
|
|
46 |
|
|
47 This work was supported by the National Science Foundation, under |
|
|
48 grants DMS-9504974 and DMS-9803599. |
|
|
49 |
|
|
50 Copyright and License: |
|
|
51 |
|
|
52 Copyright (c) 1998-2003 by the University of Florida. |
|
|
53 All Rights Reserved. |
|
|
54 |
|
|
55 THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY |
|
|
56 EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. |
|
|
57 |
|
|
58 Permission is hereby granted to use, copy, modify, and/or distribute |
|
|
59 this program, provided that the Copyright, this License, and the |
|
|
60 Availability of the original version is retained on all copies and made |
|
|
61 accessible to the end-user of any code or package that includes COLAMD |
|
|
62 or any modified version of COLAMD. |
|
|
63 |
|
|
64 Availability: |
|
|
65 |
|
|
66 The colamd/symamd library is available at |
|
|
67 |
|
|
68 http://www.cise.ufl.edu/research/sparse/colamd/ |
|
|
69 |
|
|
70 This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.c |
|
|
71 file. It requires the colamd.h file. It is required by the colamdmex.c |
|
|
72 and symamdmex.c files, for the MATLAB interface to colamd and symamd. |
|
|
73 |
|
|
74 See the ChangeLog file for changes since Version 1.0. |
|
|
75 |
|
|
76 */ |
|
|
77 |
|
|
78 /* ========================================================================== */ |
|
|
79 /* === Description of user-callable routines ================================ */ |
|
|
80 /* ========================================================================== */ |
|
|
81 |
|
|
82 /* |
|
|
83 ---------------------------------------------------------------------------- |
|
|
84 colamd_recommended: |
|
|
85 ---------------------------------------------------------------------------- |
|
|
86 |
|
|
87 C syntax: |
|
|
88 |
|
|
89 #include "colamd.h" |
|
|
90 int colamd_recommended (int nnz, int n_row, int n_col) ; |
|
|
91 |
|
|
92 or as a C macro |
|
|
93 |
|
|
94 #include "colamd.h" |
|
|
95 Alen = COLAMD_RECOMMENDED (int nnz, int n_row, int n_col) ; |
|
|
96 |
|
|
97 Purpose: |
|
|
98 |
|
|
99 Returns recommended value of Alen for use by colamd. Returns -1 |
|
|
100 if any input argument is negative. The use of this routine |
|
|
101 or macro is optional. Note that the macro uses its arguments |
|
|
102 more than once, so be careful for side effects, if you pass |
|
|
103 expressions as arguments to COLAMD_RECOMMENDED. Not needed for |
|
|
104 symamd, which dynamically allocates its own memory. |
|
|
105 |
|
|
106 Arguments (all input arguments): |
|
|
107 |
|
|
108 int nnz ; Number of nonzeros in the matrix A. This must |
|
|
109 be the same value as p [n_col] in the call to |
|
|
110 colamd - otherwise you will get a wrong value |
|
|
111 of the recommended memory to use. |
|
|
112 |
|
|
113 int n_row ; Number of rows in the matrix A. |
|
|
114 |
|
|
115 int n_col ; Number of columns in the matrix A. |
|
|
116 |
|
|
117 ---------------------------------------------------------------------------- |
|
|
118 colamd_set_defaults: |
|
|
119 ---------------------------------------------------------------------------- |
|
|
120 |
|
|
121 C syntax: |
|
|
122 |
|
|
123 #include "colamd.h" |
|
|
124 colamd_set_defaults (double knobs [COLAMD_KNOBS]) ; |
|
|
125 |
|
|
126 Purpose: |
|
|
127 |
|
|
128 Sets the default parameters. The use of this routine is optional. |
|
|
129 |
|
|
130 Arguments: |
|
|
131 |
|
|
132 double knobs [COLAMD_KNOBS] ; Output only. |
|
|
133 |
|
|
134 Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) |
|
|
135 entries are removed prior to ordering. Columns with more than |
|
|
136 (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to |
|
|
137 ordering, and placed last in the output column ordering. |
|
|
138 |
|
|
139 Symamd: uses only knobs [COLAMD_DENSE_ROW], which is knobs [0]. |
|
|
140 Rows and columns with more than (knobs [COLAMD_DENSE_ROW] * n) |
|
|
141 entries are removed prior to ordering, and placed last in the |
|
|
142 output ordering. |
|
|
143 |
|
|
144 COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1, |
|
|
145 respectively, in colamd.h. Default values of these two knobs |
|
|
146 are both 0.5. Currently, only knobs [0] and knobs [1] are |
|
|
147 used, but future versions may use more knobs. If so, they will |
|
|
148 be properly set to their defaults by the future version of |
|
|
149 colamd_set_defaults, so that the code that calls colamd will |
|
|
150 not need to change, assuming that you either use |
|
|
151 colamd_set_defaults, or pass a (double *) NULL pointer as the |
|
|
152 knobs array to colamd or symamd. |
|
|
153 |
|
|
154 ---------------------------------------------------------------------------- |
|
|
155 colamd: |
|
|
156 ---------------------------------------------------------------------------- |
|
|
157 |
|
|
158 C syntax: |
|
|
159 |
|
|
160 #include "colamd.h" |
|
|
161 int colamd (int n_row, int n_col, int Alen, int *A, int *p, |
|
|
162 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS]) ; |
|
|
163 |
|
|
164 Purpose: |
|
|
165 |
|
|
166 Computes a column ordering (Q) of A such that P(AQ)=LU or |
|
|
167 (AQ)'AQ=LL' have less fill-in and require fewer floating point |
|
|
168 operations than factorizing the unpermuted matrix A or A'A, |
|
|
169 respectively. |
|
|
170 |
|
|
171 Returns: |
|
|
172 |
|
|
173 TRUE (1) if successful, FALSE (0) otherwise. |
|
|
174 |
|
|
175 Arguments: |
|
|
176 |
|
|
177 int n_row ; Input argument. |
|
|
178 |
|
|
179 Number of rows in the matrix A. |
|
|
180 Restriction: n_row >= 0. |
|
|
181 Colamd returns FALSE if n_row is negative. |
|
|
182 |
|
|
183 int n_col ; Input argument. |
|
|
184 |
|
|
185 Number of columns in the matrix A. |
|
|
186 Restriction: n_col >= 0. |
|
|
187 Colamd returns FALSE if n_col is negative. |
|
|
188 |
|
|
189 int Alen ; Input argument. |
|
|
190 |
|
|
191 Restriction (see note): |
|
|
192 Alen >= 2*nnz + 6*(n_col+1) + 4*(n_row+1) + n_col |
|
|
193 Colamd returns FALSE if these conditions are not met. |
|
|
194 |
|
|
195 Note: this restriction makes an modest assumption regarding |
|
|
196 the size of the two typedef's structures in colamd.h. |
|
|
197 We do, however, guarantee that |
|
|
198 |
|
|
199 Alen >= colamd_recommended (nnz, n_row, n_col) |
|
|
200 |
|
|
201 or equivalently as a C preprocessor macro: |
|
|
202 |
|
|
203 Alen >= COLAMD_RECOMMENDED (nnz, n_row, n_col) |
|
|
204 |
|
|
205 will be sufficient. |
|
|
206 |
|
|
207 int A [Alen] ; Input argument, undefined on output. |
|
|
208 |
|
|
209 A is an integer array of size Alen. Alen must be at least as |
|
|
210 large as the bare minimum value given above, but this is very |
|
|
211 low, and can result in excessive run time. For best |
|
|
212 performance, we recommend that Alen be greater than or equal to |
|
|
213 colamd_recommended (nnz, n_row, n_col), which adds |
|
|
214 nnz/5 to the bare minimum value given above. |
|
|
215 |
|
|
216 On input, the row indices of the entries in column c of the |
|
|
217 matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices |
|
|
218 in a given column c need not be in ascending order, and |
|
|
219 duplicate row indices may be be present. However, colamd will |
|
|
220 work a little faster if both of these conditions are met |
|
|
221 (Colamd puts the matrix into this format, if it finds that the |
|
|
222 the conditions are not met). |
|
|
223 |
|
|
224 The matrix is 0-based. That is, rows are in the range 0 to |
|
|
225 n_row-1, and columns are in the range 0 to n_col-1. Colamd |
|
|
226 returns FALSE if any row index is out of range. |
|
|
227 |
|
|
228 The contents of A are modified during ordering, and are |
|
|
229 undefined on output. |
|
|
230 |
|
|
231 int p [n_col+1] ; Both input and output argument. |
|
|
232 |
|
|
233 p is an integer array of size n_col+1. On input, it holds the |
|
|
234 "pointers" for the column form of the matrix A. Column c of |
|
|
235 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first |
|
|
236 entry, p [0], must be zero, and p [c] <= p [c+1] must hold |
|
|
237 for all c in the range 0 to n_col-1. The value p [n_col] is |
|
|
238 thus the total number of entries in the pattern of the matrix A. |
|
|
239 Colamd returns FALSE if these conditions are not met. |
|
|
240 |
|
|
241 On output, if colamd returns TRUE, the array p holds the column |
|
|
242 permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is |
|
|
243 the first column index in the new ordering, and p [n_col-1] is |
|
|
244 the last. That is, p [k] = j means that column j of A is the |
|
|
245 kth pivot column, in AQ, where k is in the range 0 to n_col-1 |
|
|
246 (p [0] = j means that column j of A is the first column in AQ). |
|
|
247 |
|
|
248 If colamd returns FALSE, then no permutation is returned, and |
|
|
249 p is undefined on output. |
|
|
250 |
|
|
251 double knobs [COLAMD_KNOBS] ; Input argument. |
|
|
252 |
|
|
253 See colamd_set_defaults for a description. |
|
|
254 |
|
|
255 int stats [COLAMD_STATS] ; Output argument. |
|
|
256 |
|
|
257 Statistics on the ordering, and error status. |
|
|
258 See colamd.h for related definitions. |
|
|
259 Colamd returns FALSE if stats is not present. |
|
|
260 |
|
|
261 stats [0]: number of dense or empty rows ignored. |
|
|
262 |
|
|
263 stats [1]: number of dense or empty columns ignored (and |
|
|
264 ordered last in the output permutation p) |
|
|
265 Note that a row can become "empty" if it |
|
|
266 contains only "dense" and/or "empty" columns, |
|
|
267 and similarly a column can become "empty" if it |
|
|
268 only contains "dense" and/or "empty" rows. |
|
|
269 |
|
|
270 stats [2]: number of garbage collections performed. |
|
|
271 This can be excessively high if Alen is close |
|
|
272 to the minimum required value. |
|
|
273 |
|
|
274 stats [3]: status code. < 0 is an error code. |
|
|
275 > 1 is a warning or notice. |
|
|
276 |
|
|
277 0 OK. Each column of the input matrix contained |
|
|
278 row indices in increasing order, with no |
|
|
279 duplicates. |
|
|
280 |
|
|
281 1 OK, but columns of input matrix were jumbled |
|
|
282 (unsorted columns or duplicate entries). Colamd |
|
|
283 had to do some extra work to sort the matrix |
|
|
284 first and remove duplicate entries, but it |
|
|
285 still was able to return a valid permutation |
|
|
286 (return value of colamd was TRUE). |
|
|
287 |
|
|
288 stats [4]: highest numbered column that |
|
|
289 is unsorted or has duplicate |
|
|
290 entries. |
|
|
291 stats [5]: last seen duplicate or |
|
|
292 unsorted row index. |
|
|
293 stats [6]: number of duplicate or |
|
|
294 unsorted row indices. |
|
|
295 |
|
|
296 -1 A is a null pointer |
|
|
297 |
|
|
298 -2 p is a null pointer |
|
|
299 |
|
|
300 -3 n_row is negative |
|
|
301 |
|
|
302 stats [4]: n_row |
|
|
303 |
|
|
304 -4 n_col is negative |
|
|
305 |
|
|
306 stats [4]: n_col |
|
|
307 |
|
|
308 -5 number of nonzeros in matrix is negative |
|
|
309 |
|
|
310 stats [4]: number of nonzeros, p [n_col] |
|
|
311 |
|
|
312 -6 p [0] is nonzero |
|
|
313 |
|
|
314 stats [4]: p [0] |
|
|
315 |
|
|
316 -7 A is too small |
|
|
317 |
|
|
318 stats [4]: required size |
|
|
319 stats [5]: actual size (Alen) |
|
|
320 |
|
|
321 -8 a column has a negative number of entries |
|
|
322 |
|
|
323 stats [4]: column with < 0 entries |
|
|
324 stats [5]: number of entries in col |
|
|
325 |
|
|
326 -9 a row index is out of bounds |
|
|
327 |
|
|
328 stats [4]: column with bad row index |
|
|
329 stats [5]: bad row index |
|
|
330 stats [6]: n_row, # of rows of matrx |
|
|
331 |
|
|
332 -10 (unused; see symamd.c) |
|
|
333 |
|
|
334 -999 (unused; see symamd.c) |
|
|
335 |
|
|
336 Future versions may return more statistics in the stats array. |
|
|
337 |
|
|
338 Example: |
|
|
339 |
|
|
340 See http://www.cise.ufl.edu/research/sparse/colamd/example.c |
|
|
341 for a complete example. |
|
|
342 |
|
|
343 To order the columns of a 5-by-4 matrix with 11 nonzero entries in |
|
|
344 the following nonzero pattern |
|
|
345 |
|
|
346 x 0 x 0 |
|
|
347 x 0 x x |
|
|
348 0 x x 0 |
|
|
349 0 0 x x |
|
|
350 x x 0 0 |
|
|
351 |
|
|
352 with default knobs and no output statistics, do the following: |
|
|
353 |
|
|
354 #include "colamd.h" |
|
|
355 #define ALEN COLAMD_RECOMMENDED (11, 5, 4) |
|
|
356 int A [ALEN] = {1, 2, 5, 3, 5, 1, 2, 3, 4, 2, 4} ; |
|
|
357 int p [ ] = {0, 3, 5, 9, 11} ; |
|
|
358 int stats [COLAMD_STATS] ; |
|
|
359 colamd (5, 4, ALEN, A, p, (double *) NULL, stats) ; |
|
|
360 |
|
|
361 The permutation is returned in the array p, and A is destroyed. |
|
|
362 |
|
|
363 ---------------------------------------------------------------------------- |
|
|
364 symamd: |
|
|
365 ---------------------------------------------------------------------------- |
|
|
366 |
|
|
367 C syntax: |
|
|
368 |
|
|
369 #include "colamd.h" |
|
|
370 int symamd (int n, int *A, int *p, int *perm, |
|
|
371 double knobs [COLAMD_KNOBS], int stats [COLAMD_STATS], |
|
|
372 void (*allocate) (size_t, size_t), void (*release) (void *)) ; |
|
|
373 |
|
|
374 Purpose: |
|
|
375 |
|
|
376 The symamd routine computes an ordering P of a symmetric sparse |
|
|
377 matrix A such that the Cholesky factorization PAP' = LL' remains |
|
|
378 sparse. It is based on a column ordering of a matrix M constructed |
|
|
379 so that the nonzero pattern of M'M is the same as A. The matrix A |
|
|
380 is assumed to be symmetric; only the strictly lower triangular part |
|
|
381 is accessed. You must pass your selected memory allocator (usually |
|
|
382 calloc/free or mxCalloc/mxFree) to symamd, for it to allocate |
|
|
383 memory for the temporary matrix M. |
|
|
384 |
|
|
385 Returns: |
|
|
386 |
|
|
387 TRUE (1) if successful, FALSE (0) otherwise. |
|
|
388 |
|
|
389 Arguments: |
|
|
390 |
|
|
391 int n ; Input argument. |
|
|
392 |
|
|
393 Number of rows and columns in the symmetrix matrix A. |
|
|
394 Restriction: n >= 0. |
|
|
395 Symamd returns FALSE if n is negative. |
|
|
396 |
|
|
397 int A [nnz] ; Input argument. |
|
|
398 |
|
|
399 A is an integer array of size nnz, where nnz = p [n]. |
|
|
400 |
|
|
401 The row indices of the entries in column c of the matrix are |
|
|
402 held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a |
|
|
403 given column c need not be in ascending order, and duplicate |
|
|
404 row indices may be present. However, symamd will run faster |
|
|
405 if the columns are in sorted order with no duplicate entries. |
|
|
406 |
|
|
407 The matrix is 0-based. That is, rows are in the range 0 to |
|
|
408 n-1, and columns are in the range 0 to n-1. Symamd |
|
|
409 returns FALSE if any row index is out of range. |
|
|
410 |
|
|
411 The contents of A are not modified. |
|
|
412 |
|
|
413 int p [n+1] ; Input argument. |
|
|
414 |
|
|
415 p is an integer array of size n+1. On input, it holds the |
|
|
416 "pointers" for the column form of the matrix A. Column c of |
|
|
417 the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first |
|
|
418 entry, p [0], must be zero, and p [c] <= p [c+1] must hold |
|
|
419 for all c in the range 0 to n-1. The value p [n] is |
|
|
420 thus the total number of entries in the pattern of the matrix A. |
|
|
421 Symamd returns FALSE if these conditions are not met. |
|
|
422 |
|
|
423 The contents of p are not modified. |
|
|
424 |
|
|
425 int perm [n+1] ; Output argument. |
|
|
426 |
|
|
427 On output, if symamd returns TRUE, the array perm holds the |
|
|
428 permutation P, where perm [0] is the first index in the new |
|
|
429 ordering, and perm [n-1] is the last. That is, perm [k] = j |
|
|
430 means that row and column j of A is the kth column in PAP', |
|
|
431 where k is in the range 0 to n-1 (perm [0] = j means |
|
|
432 that row and column j of A are the first row and column in |
|
|
433 PAP'). The array is used as a workspace during the ordering, |
|
|
434 which is why it must be of length n+1, not just n. |
|
|
435 |
|
|
436 double knobs [COLAMD_KNOBS] ; Input argument. |
|
|
437 |
|
|
438 See colamd_set_defaults for a description. |
|
|
439 |
|
|
440 int stats [COLAMD_STATS] ; Output argument. |
|
|
441 |
|
|
442 Statistics on the ordering, and error status. |
|
|
443 See colamd.h for related definitions. |
|
|
444 Symamd returns FALSE if stats is not present. |
|
|
445 |
|
|
446 stats [0]: number of dense or empty row and columns ignored |
|
|
447 (and ordered last in the output permutation |
|
|
448 perm). Note that a row/column can become |
|
|
449 "empty" if it contains only "dense" and/or |
|
|
450 "empty" columns/rows. |
|
|
451 |
|
|
452 stats [1]: (same as stats [0]) |
|
|
453 |
|
|
454 stats [2]: number of garbage collections performed. |
|
|
455 |
|
|
456 stats [3]: status code. < 0 is an error code. |
|
|
457 > 1 is a warning or notice. |
|
|
458 |
|
|
459 0 OK. Each column of the input matrix contained |
|
|
460 row indices in increasing order, with no |
|
|
461 duplicates. |
|
|
462 |
|
|
463 1 OK, but columns of input matrix were jumbled |
|
|
464 (unsorted columns or duplicate entries). Symamd |
|
|
465 had to do some extra work to sort the matrix |
|
|
466 first and remove duplicate entries, but it |
|
|
467 still was able to return a valid permutation |
|
|
468 (return value of symamd was TRUE). |
|
|
469 |
|
|
470 stats [4]: highest numbered column that |
|
|
471 is unsorted or has duplicate |
|
|
472 entries. |
|
|
473 stats [5]: last seen duplicate or |
|
|
474 unsorted row index. |
|
|
475 stats [6]: number of duplicate or |
|
|
476 unsorted row indices. |
|
|
477 |
|
|
478 -1 A is a null pointer |
|
|
479 |
|
|
480 -2 p is a null pointer |
|
|
481 |
|
|
482 -3 (unused, see colamd.c) |
|
|
483 |
|
|
484 -4 n is negative |
|
|
485 |
|
|
486 stats [4]: n |
|
|
487 |
|
|
488 -5 number of nonzeros in matrix is negative |
|
|
489 |
|
|
490 stats [4]: # of nonzeros (p [n]). |
|
|
491 |
|
|
492 -6 p [0] is nonzero |
|
|
493 |
|
|
494 stats [4]: p [0] |
|
|
495 |
|
|
496 -7 (unused) |
|
|
497 |
|
|
498 -8 a column has a negative number of entries |
|
|
499 |
|
|
500 stats [4]: column with < 0 entries |
|
|
501 stats [5]: number of entries in col |
|
|
502 |
|
|
503 -9 a row index is out of bounds |
|
|
504 |
|
|
505 stats [4]: column with bad row index |
|
|
506 stats [5]: bad row index |
|
|
507 stats [6]: n_row, # of rows of matrx |
|
|
508 |
|
|
509 -10 out of memory (unable to allocate temporary |
|
|
510 workspace for M or count arrays using the |
|
|
511 "allocate" routine passed into symamd). |
|
|
512 |
|
|
513 -999 internal error. colamd failed to order the |
|
|
514 matrix M, when it should have succeeded. This |
|
|
515 indicates a bug. If this (and *only* this) |
|
|
516 error code occurs, please contact the authors. |
|
|
517 Don't contact the authors if you get any other |
|
|
518 error code. |
|
|
519 |
|
|
520 Future versions may return more statistics in the stats array. |
|
|
521 |
|
|
522 void * (*allocate) (size_t, size_t) |
|
|
523 |
|
|
524 A pointer to a function providing memory allocation. The |
|
|
525 allocated memory must be returned initialized to zero. For a |
|
|
526 C application, this argument should normally be a pointer to |
|
|
527 calloc. For a MATLAB mexFunction, the routine mxCalloc is |
|
|
528 passed instead. |
|
|
529 |
|
|
530 void (*release) (size_t, size_t) |
|
|
531 |
|
|
532 A pointer to a function that frees memory allocated by the |
|
|
533 memory allocation routine above. For a C application, this |
|
|
534 argument should normally be a pointer to free. For a MATLAB |
|
|
535 mexFunction, the routine mxFree is passed instead. |
|
|
536 |
|
|
537 |
|
|
538 ---------------------------------------------------------------------------- |
|
|
539 colamd_report: |
|
|
540 ---------------------------------------------------------------------------- |
|
|
541 |
|
|
542 C syntax: |
|
|
543 |
|
|
544 #include "colamd.h" |
|
|
545 colamd_report (int stats [COLAMD_STATS]) ; |
|
|
546 |
|
|
547 Purpose: |
|
|
548 |
|
|
549 Prints the error status and statistics recorded in the stats |
|
|
550 array on the standard error output (for a standard C routine) |
|
|
551 or on the MATLAB output (for a mexFunction). |
|
|
552 |
|
|
553 Arguments: |
|
|
554 |
|
|
555 int stats [COLAMD_STATS] ; Input only. Statistics from colamd. |
|
|
556 |
|
|
557 |
|
|
558 ---------------------------------------------------------------------------- |
|
|
559 symamd_report: |
|
|
560 ---------------------------------------------------------------------------- |
|
|
561 |
|
|
562 C syntax: |
|
|
563 |
|
|
564 #include "colamd.h" |
|
|
565 symamd_report (int stats [COLAMD_STATS]) ; |
|
|
566 |
|
|
567 Purpose: |
|
|
568 |
|
|
569 Prints the error status and statistics recorded in the stats |
|
|
570 array on the standard error output (for a standard C routine) |
|
|
571 or on the MATLAB output (for a mexFunction). |
|
|
572 |
|
|
573 Arguments: |
|
|
574 |
|
|
575 int stats [COLAMD_STATS] ; Input only. Statistics from symamd. |
|
|
576 |
|
|
577 |
|
|
578 */ |
|
|
579 |
|
|
580 /* ========================================================================== */ |
|
|
581 /* === Scaffolding code definitions ======================================== */ |
|
|
582 /* ========================================================================== */ |
|
|
583 |
|
|
584 /* Ensure that debugging is turned off: */ |
|
|
585 #ifndef NDEBUG |
|
|
586 #define NDEBUG |
|
|
587 #endif /* NDEBUG */ |
|
|
588 |
|
|
589 /* |
|
|
590 Our "scaffolding code" philosophy: In our opinion, well-written library |
|
|
591 code should keep its "debugging" code, and just normally have it turned off |
|
|
592 by the compiler so as not to interfere with performance. This serves |
|
|
593 several purposes: |
|
|
594 |
|
|
595 (1) assertions act as comments to the reader, telling you what the code |
|
|
596 expects at that point. All assertions will always be true (unless |
|
|
597 there really is a bug, of course). |
|
|
598 |
|
|
599 (2) leaving in the scaffolding code assists anyone who would like to modify |
|
|
600 the code, or understand the algorithm (by reading the debugging output, |
|
|
601 one can get a glimpse into what the code is doing). |
|
|
602 |
|
|
603 (3) (gasp!) for actually finding bugs. This code has been heavily tested |
|
|
604 and "should" be fully functional and bug-free ... but you never know... |
|
|
605 |
|
|
606 To enable debugging, comment out the "#define NDEBUG" above. For a MATLAB |
|
|
607 mexFunction, you will also need to modify mexopts.sh to remove the -DNDEBUG |
|
|
608 definition. The code will become outrageously slow when debugging is |
|
|
609 enabled. To control the level of debugging output, set an environment |
|
|
610 variable D to 0 (little), 1 (some), 2, 3, or 4 (lots). When debugging, |
|
|
611 you should see the following message on the standard output: |
|
|
612 |
|
|
613 colamd: debug version, D = 1 (THIS WILL BE SLOW!) |
|
|
614 |
|
|
615 or a similar message for symamd. If you don't, then debugging has not |
|
|
616 been enabled. |
|
|
617 |
|
|
618 */ |
|
|
619 |
|
|
620 /* ========================================================================== */ |
|
|
621 /* === Include files ======================================================== */ |
|
|
622 /* ========================================================================== */ |
|
|
623 |
|
|
624 #include "colamd.h" |
|
|
625 #include <limits.h> |
|
|
626 |
|
|
627 #ifdef MATLAB_MEX_FILE |
|
|
628 #include "mex.h" |
|
|
629 #include "matrix.h" |
|
|
630 #else |
|
|
631 #include <stdio.h> |
|
|
632 #include <assert.h> |
|
|
633 #endif /* MATLAB_MEX_FILE */ |
|
|
634 |
|
|
635 /* ========================================================================== */ |
|
|
636 /* === Definitions ========================================================== */ |
|
|
637 /* ========================================================================== */ |
|
|
638 |
|
|
639 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */ |
|
|
640 #define PUBLIC |
|
|
641 #define PRIVATE static |
|
|
642 |
|
|
643 #define MAX(a,b) (((a) > (b)) ? (a) : (b)) |
|
|
644 #define MIN(a,b) (((a) < (b)) ? (a) : (b)) |
|
|
645 |
|
|
646 #define ONES_COMPLEMENT(r) (-(r)-1) |
|
|
647 |
|
|
648 /* -------------------------------------------------------------------------- */ |
|
|
649 /* Change for version 2.1: define TRUE and FALSE only if not yet defined */ |
|
|
650 /* -------------------------------------------------------------------------- */ |
|
|
651 |
|
|
652 #ifndef TRUE |
|
|
653 #define TRUE (1) |
|
|
654 #endif |
|
|
655 |
|
|
656 #ifndef FALSE |
|
|
657 #define FALSE (0) |
|
|
658 #endif |
|
|
659 |
|
|
660 /* -------------------------------------------------------------------------- */ |
|
|
661 |
|
|
662 #define EMPTY (-1) |
|
|
663 |
|
|
664 /* Row and column status */ |
|
|
665 #define ALIVE (0) |
|
|
666 #define DEAD (-1) |
|
|
667 |
|
|
668 /* Column status */ |
|
|
669 #define DEAD_PRINCIPAL (-1) |
|
|
670 #define DEAD_NON_PRINCIPAL (-2) |
|
|
671 |
|
|
672 /* Macros for row and column status update and checking. */ |
|
|
673 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) |
|
|
674 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) |
|
|
675 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) |
|
|
676 #define COL_IS_DEAD(c) (Col [c].start < ALIVE) |
|
|
677 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) |
|
|
678 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) |
|
|
679 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } |
|
|
680 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } |
|
|
681 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } |
|
|
682 |
|
|
683 /* ========================================================================== */ |
|
|
684 /* === Colamd reporting mechanism =========================================== */ |
|
|
685 /* ========================================================================== */ |
|
|
686 |
|
|
687 #ifdef MATLAB_MEX_FILE |
|
|
688 |
|
|
689 /* use mexPrintf in a MATLAB mexFunction, for debugging and statistics output */ |
|
|
690 #define PRINTF mexPrintf |
|
|
691 |
|
|
692 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */ |
|
|
693 #define INDEX(i) ((i)+1) |
|
|
694 |
|
|
695 #else |
|
|
696 |
|
|
697 /* Use printf in standard C environment, for debugging and statistics output. */ |
|
|
698 /* Output is generated only if debugging is enabled at compile time, or if */ |
|
|
699 /* the caller explicitly calls colamd_report or symamd_report. */ |
|
|
700 #define PRINTF printf |
|
|
701 |
|
|
702 /* In C, matrices are 0-based and indices are reported as such in *_report */ |
|
|
703 #define INDEX(i) (i) |
|
|
704 |
|
|
705 #endif /* MATLAB_MEX_FILE */ |
|
|
706 |
|
|
707 /* ========================================================================== */ |
|
|
708 /* === Prototypes of PRIVATE routines ======================================= */ |
|
|
709 /* ========================================================================== */ |
|
|
710 |
|
|
711 PRIVATE int init_rows_cols |
|
|
712 ( |
|
|
713 int n_row, |
|
|
714 int n_col, |
|
|
715 Colamd_Row Row [], |
|
|
716 Colamd_Col Col [], |
|
|
717 int A [], |
|
|
718 int p [], |
|
|
719 int stats [COLAMD_STATS] |
|
|
720 ) ; |
|
|
721 |
|
|
722 PRIVATE void init_scoring |
|
|
723 ( |
|
|
724 int n_row, |
|
|
725 int n_col, |
|
|
726 Colamd_Row Row [], |
|
|
727 Colamd_Col Col [], |
|
|
728 int A [], |
|
|
729 int head [], |
|
|
730 double knobs [COLAMD_KNOBS], |
|
|
731 int *p_n_row2, |
|
|
732 int *p_n_col2, |
|
|
733 int *p_max_deg |
|
|
734 ) ; |
|
|
735 |
|
|
736 PRIVATE int find_ordering |
|
|
737 ( |
|
|
738 int n_row, |
|
|
739 int n_col, |
|
|
740 int Alen, |
|
|
741 Colamd_Row Row [], |
|
|
742 Colamd_Col Col [], |
|
|
743 int A [], |
|
|
744 int head [], |
|
|
745 int n_col2, |
|
|
746 int max_deg, |
|
|
747 int pfree |
|
|
748 ) ; |
|
|
749 |
|
|
750 PRIVATE void order_children |
|
|
751 ( |
|
|
752 int n_col, |
|
|
753 Colamd_Col Col [], |
|
|
754 int p [] |
|
|
755 ) ; |
|
|
756 |
|
|
757 PRIVATE void detect_super_cols |
|
|
758 ( |
|
|
759 |
|
|
760 #ifndef NDEBUG |
|
|
761 int n_col, |
|
|
762 Colamd_Row Row [], |
|
|
763 #endif /* NDEBUG */ |
|
|
764 |
|
|
765 Colamd_Col Col [], |
|
|
766 int A [], |
|
|
767 int head [], |
|
|
768 int row_start, |
|
|
769 int row_length |
|
|
770 ) ; |
|
|
771 |
|
|
772 PRIVATE int garbage_collection |
|
|
773 ( |
|
|
774 int n_row, |
|
|
775 int n_col, |
|
|
776 Colamd_Row Row [], |
|
|
777 Colamd_Col Col [], |
|
|
778 int A [], |
|
|
779 int *pfree |
|
|
780 ) ; |
|
|
781 |
|
|
782 PRIVATE int clear_mark |
|
|
783 ( |
|
|
784 int n_row, |
|
|
785 Colamd_Row Row [] |
|
|
786 ) ; |
|
|
787 |
|
|
788 PRIVATE void print_report |
|
|
789 ( |
|
|
790 char *method, |
|
|
791 int stats [COLAMD_STATS] |
|
|
792 ) ; |
|
|
793 |
|
|
794 /* ========================================================================== */ |
|
|
795 /* === Debugging prototypes and definitions ================================= */ |
|
|
796 /* ========================================================================== */ |
|
|
797 |
|
|
798 #ifndef NDEBUG |
|
|
799 |
|
|
800 /* colamd_debug is the *ONLY* global variable, and is only */ |
|
|
801 /* present when debugging */ |
|
|
802 |
|
|
803 PRIVATE int colamd_debug ; /* debug print level */ |
|
|
804 |
|
|
805 #define DEBUG0(params) { (void) PRINTF params ; } |
|
|
806 #define DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; } |
|
|
807 #define DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; } |
|
|
808 #define DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; } |
|
|
809 #define DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; } |
|
|
810 |
|
|
811 #ifdef MATLAB_MEX_FILE |
|
|
812 #define ASSERT(expression) (mxAssert ((expression), "")) |
|
|
813 #else |
|
|
814 #define ASSERT(expression) (assert (expression)) |
|
|
815 #endif /* MATLAB_MEX_FILE */ |
|
|
816 |
|
|
817 PRIVATE void colamd_get_debug /* gets the debug print level from getenv */ |
|
|
818 ( |
|
|
819 char *method |
|
|
820 ) ; |
|
|
821 |
|
|
822 PRIVATE void debug_deg_lists |
|
|
823 ( |
|
|
824 int n_row, |
|
|
825 int n_col, |
|
|
826 Colamd_Row Row [], |
|
|
827 Colamd_Col Col [], |
|
|
828 int head [], |
|
|
829 int min_score, |
|
|
830 int should, |
|
|
831 int max_deg |
|
|
832 ) ; |
|
|
833 |
|
|
834 PRIVATE void debug_mark |
|
|
835 ( |
|
|
836 int n_row, |
|
|
837 Colamd_Row Row [], |
|
|
838 int tag_mark, |
|
|
839 int max_mark |
|
|
840 ) ; |
|
|
841 |
|
|
842 PRIVATE void debug_matrix |
|
|
843 ( |
|
|
844 int n_row, |
|
|
845 int n_col, |
|
|
846 Colamd_Row Row [], |
|
|
847 Colamd_Col Col [], |
|
|
848 int A [] |
|
|
849 ) ; |
|
|
850 |
|
|
851 PRIVATE void debug_structures |
|
|
852 ( |
|
|
853 int n_row, |
|
|
854 int n_col, |
|
|
855 Colamd_Row Row [], |
|
|
856 Colamd_Col Col [], |
|
|
857 int A [], |
|
|
858 int n_col2 |
|
|
859 ) ; |
|
|
860 |
|
|
861 #else /* NDEBUG */ |
|
|
862 |
|
|
863 /* === No debugging ========================================================= */ |
|
|
864 |
|
|
865 #define DEBUG0(params) ; |
|
|
866 #define DEBUG1(params) ; |
|
|
867 #define DEBUG2(params) ; |
|
|
868 #define DEBUG3(params) ; |
|
|
869 #define DEBUG4(params) ; |
|
|
870 |
|
|
871 #define ASSERT(expression) ((void) 0) |
|
|
872 |
|
|
873 #endif /* NDEBUG */ |
|
|
874 |
|
|
875 /* ========================================================================== */ |
|
|
876 |
|
|
877 |
|
|
878 |
|
|
879 /* ========================================================================== */ |
|
|
880 /* === USER-CALLABLE ROUTINES: ============================================== */ |
|
|
881 /* ========================================================================== */ |
|
|
882 |
|
|
883 |
|
|
884 /* ========================================================================== */ |
|
|
885 /* === colamd_recommended =================================================== */ |
|
|
886 /* ========================================================================== */ |
|
|
887 |
|
|
888 /* |
|
|
889 The colamd_recommended routine returns the suggested size for Alen. This |
|
|
890 value has been determined to provide good balance between the number of |
|
|
891 garbage collections and the memory requirements for colamd. If any |
|
|
892 argument is negative, a -1 is returned as an error condition. This |
|
|
893 function is also available as a macro defined in colamd.h, so that you |
|
|
894 can use it for a statically-allocated array size. |
|
|
895 */ |
|
|
896 |
|
|
897 PUBLIC int colamd_recommended /* returns recommended value of Alen. */ |
|
|
898 ( |
|
|
899 /* === Parameters ======================================================= */ |
|
|
900 |
|
|
901 int nnz, /* number of nonzeros in A */ |
|
|
902 int n_row, /* number of rows in A */ |
|
|
903 int n_col /* number of columns in A */ |
|
|
904 ) |
|
|
905 { |
|
|
906 return (COLAMD_RECOMMENDED (nnz, n_row, n_col)) ; |
|
|
907 } |
|
|
908 |
|
|
909 |
|
|
910 /* ========================================================================== */ |
|
|
911 /* === colamd_set_defaults ================================================== */ |
|
|
912 /* ========================================================================== */ |
|
|
913 |
|
|
914 /* |
|
|
915 The colamd_set_defaults routine sets the default values of the user- |
|
|
916 controllable parameters for colamd: |
|
|
917 |
|
|
918 knobs [0] rows with knobs[0]*n_col entries or more are removed |
|
|
919 prior to ordering in colamd. Rows and columns with |
|
|
920 knobs[0]*n_col entries or more are removed prior to |
|
|
921 ordering in symamd and placed last in the output |
|
|
922 ordering. |
|
|
923 |
|
|
924 knobs [1] columns with knobs[1]*n_row entries or more are removed |
|
|
925 prior to ordering in colamd, and placed last in the |
|
|
926 column permutation. Symamd ignores this knob. |
|
|
927 |
|
|
928 knobs [2..19] unused, but future versions might use this |
|
|
929 */ |
|
|
930 |
|
|
931 PUBLIC void colamd_set_defaults |
|
|
932 ( |
|
|
933 /* === Parameters ======================================================= */ |
|
|
934 |
|
|
935 double knobs [COLAMD_KNOBS] /* knob array */ |
|
|
936 ) |
|
|
937 { |
|
|
938 /* === Local variables ================================================== */ |
|
|
939 |
|
|
940 int i ; |
|
|
941 |
|
|
942 if (!knobs) |
|
|
943 { |
|
|
944 return ; /* no knobs to initialize */ |
|
|
945 } |
|
|
946 for (i = 0 ; i < COLAMD_KNOBS ; i++) |
|
|
947 { |
|
|
948 knobs [i] = 0 ; |
|
|
949 } |
|
|
950 knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ |
|
|
951 knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ |
|
|
952 } |
|
|
953 |
|
|
954 |
|
|
955 /* ========================================================================== */ |
|
|
956 /* === symamd =============================================================== */ |
|
|
957 /* ========================================================================== */ |
|
|
958 |
|
|
959 PUBLIC int symamd /* return TRUE if OK, FALSE otherwise */ |
|
|
960 ( |
|
|
961 /* === Parameters ======================================================= */ |
|
|
962 |
|
|
963 int n, /* number of rows and columns of A */ |
|
|
964 int A [], /* row indices of A */ |
|
|
965 int p [], /* column pointers of A */ |
|
|
966 int perm [], /* output permutation, size n+1 */ |
|
|
967 double knobs [COLAMD_KNOBS], /* parameters (uses defaults if NULL) */ |
|
|
968 int stats [COLAMD_STATS], /* output statistics and error codes */ |
|
|
969 void * (*allocate) (size_t, size_t), |
|
|
970 /* pointer to calloc (ANSI C) or */ |
|
|
971 /* mxCalloc (for MATLAB mexFunction) */ |
|
|
972 void (*release) (void *) |
|
|
973 /* pointer to free (ANSI C) or */ |
|
|
974 /* mxFree (for MATLAB mexFunction) */ |
|
|
975 ) |
|
|
976 { |
|
|
977 /* === Local variables ================================================== */ |
|
|
978 |
|
|
979 int *count ; /* length of each column of M, and col pointer*/ |
|
|
980 int *mark ; /* mark array for finding duplicate entries */ |
|
|
981 int *M ; /* row indices of matrix M */ |
|
|
982 int Mlen ; /* length of M */ |
|
|
983 int n_row ; /* number of rows in M */ |
|
|
984 int nnz ; /* number of entries in A */ |
|
|
985 int i ; /* row index of A */ |
|
|
986 int j ; /* column index of A */ |
|
|
987 int k ; /* row index of M */ |
|
|
988 int mnz ; /* number of nonzeros in M */ |
|
|
989 int pp ; /* index into a column of A */ |
|
|
990 int last_row ; /* last row seen in the current column */ |
|
|
991 int length ; /* number of nonzeros in a column */ |
|
|
992 |
|
|
993 double cknobs [COLAMD_KNOBS] ; /* knobs for colamd */ |
|
|
994 double default_knobs [COLAMD_KNOBS] ; /* default knobs for colamd */ |
|
|
995 int cstats [COLAMD_STATS] ; /* colamd stats */ |
|
|
996 |
|
|
997 #ifndef NDEBUG |
|
|
998 colamd_get_debug ("symamd") ; |
|
|
999 #endif /* NDEBUG */ |
|
|
1000 |
|
|
1001 /* === Check the input arguments ======================================== */ |
|
|
1002 |
|
|
1003 if (!stats) |
|
|
1004 { |
|
|
1005 DEBUG0 (("symamd: stats not present\n")) ; |
|
|
1006 return (FALSE) ; |
|
|
1007 } |
|
|
1008 for (i = 0 ; i < COLAMD_STATS ; i++) |
|
|
1009 { |
|
|
1010 stats [i] = 0 ; |
|
|
1011 } |
|
|
1012 stats [COLAMD_STATUS] = COLAMD_OK ; |
|
|
1013 stats [COLAMD_INFO1] = -1 ; |
|
|
1014 stats [COLAMD_INFO2] = -1 ; |
|
|
1015 |
|
|
1016 if (!A) |
|
|
1017 { |
|
|
1018 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; |
|
|
1019 DEBUG0 (("symamd: A not present\n")) ; |
|
|
1020 return (FALSE) ; |
|
|
1021 } |
|
|
1022 |
|
|
1023 if (!p) /* p is not present */ |
|
|
1024 { |
|
|
1025 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; |
|
|
1026 DEBUG0 (("symamd: p not present\n")) ; |
|
|
1027 return (FALSE) ; |
|
|
1028 } |
|
|
1029 |
|
|
1030 if (n < 0) /* n must be >= 0 */ |
|
|
1031 { |
|
|
1032 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; |
|
|
1033 stats [COLAMD_INFO1] = n ; |
|
|
1034 DEBUG0 (("symamd: n negative %d\n", n)) ; |
|
|
1035 return (FALSE) ; |
|
|
1036 } |
|
|
1037 |
|
|
1038 nnz = p [n] ; |
|
|
1039 if (nnz < 0) /* nnz must be >= 0 */ |
|
|
1040 { |
|
|
1041 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; |
|
|
1042 stats [COLAMD_INFO1] = nnz ; |
|
|
1043 DEBUG0 (("symamd: number of entries negative %d\n", nnz)) ; |
|
|
1044 return (FALSE) ; |
|
|
1045 } |
|
|
1046 |
|
|
1047 if (p [0] != 0) |
|
|
1048 { |
|
|
1049 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; |
|
|
1050 stats [COLAMD_INFO1] = p [0] ; |
|
|
1051 DEBUG0 (("symamd: p[0] not zero %d\n", p [0])) ; |
|
|
1052 return (FALSE) ; |
|
|
1053 } |
|
|
1054 |
|
|
1055 /* === If no knobs, set default knobs =================================== */ |
|
|
1056 |
|
|
1057 if (!knobs) |
|
|
1058 { |
|
|
1059 colamd_set_defaults (default_knobs) ; |
|
|
1060 knobs = default_knobs ; |
|
|
1061 } |
|
|
1062 |
|
|
1063 /* === Allocate count and mark ========================================== */ |
|
|
1064 |
|
|
1065 count = (int *) ((*allocate) (n+1, sizeof (int))) ; |
|
|
1066 if (!count) |
|
|
1067 { |
|
|
1068 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; |
|
|
1069 DEBUG0 (("symamd: allocate count (size %d) failed\n", n+1)) ; |
|
|
1070 return (FALSE) ; |
|
|
1071 } |
|
|
1072 |
|
|
1073 mark = (int *) ((*allocate) (n+1, sizeof (int))) ; |
|
|
1074 if (!mark) |
|
|
1075 { |
|
|
1076 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; |
|
|
1077 (*release) ((void *) count) ; |
|
|
1078 DEBUG0 (("symamd: allocate mark (size %d) failed\n", n+1)) ; |
|
|
1079 return (FALSE) ; |
|
|
1080 } |
|
|
1081 |
|
|
1082 /* === Compute column counts of M, check if A is valid ================== */ |
|
|
1083 |
|
|
1084 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ |
|
|
1085 |
|
|
1086 for (i = 0 ; i < n ; i++) |
|
|
1087 { |
|
|
1088 mark [i] = -1 ; |
|
|
1089 } |
|
|
1090 |
|
|
1091 for (j = 0 ; j < n ; j++) |
|
|
1092 { |
|
|
1093 last_row = -1 ; |
|
|
1094 |
|
|
1095 length = p [j+1] - p [j] ; |
|
|
1096 if (length < 0) |
|
|
1097 { |
|
|
1098 /* column pointers must be non-decreasing */ |
|
|
1099 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; |
|
|
1100 stats [COLAMD_INFO1] = j ; |
|
|
1101 stats [COLAMD_INFO2] = length ; |
|
|
1102 (*release) ((void *) count) ; |
|
|
1103 (*release) ((void *) mark) ; |
|
|
1104 DEBUG0 (("symamd: col %d negative length %d\n", j, length)) ; |
|
|
1105 return (FALSE) ; |
|
|
1106 } |
|
|
1107 |
|
|
1108 for (pp = p [j] ; pp < p [j+1] ; pp++) |
|
|
1109 { |
|
|
1110 i = A [pp] ; |
|
|
1111 if (i < 0 || i >= n) |
|
|
1112 { |
|
|
1113 /* row index i, in column j, is out of bounds */ |
|
|
1114 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; |
|
|
1115 stats [COLAMD_INFO1] = j ; |
|
|
1116 stats [COLAMD_INFO2] = i ; |
|
|
1117 stats [COLAMD_INFO3] = n ; |
|
|
1118 (*release) ((void *) count) ; |
|
|
1119 (*release) ((void *) mark) ; |
|
|
1120 DEBUG0 (("symamd: row %d col %d out of bounds\n", i, j)) ; |
|
|
1121 return (FALSE) ; |
|
|
1122 } |
|
|
1123 |
|
|
1124 if (i <= last_row || mark [i] == j) |
|
|
1125 { |
|
|
1126 /* row index is unsorted or repeated (or both), thus col */ |
|
|
1127 /* is jumbled. This is a notice, not an error condition. */ |
|
|
1128 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; |
|
|
1129 stats [COLAMD_INFO1] = j ; |
|
|
1130 stats [COLAMD_INFO2] = i ; |
|
|
1131 (stats [COLAMD_INFO3]) ++ ; |
|
|
1132 DEBUG1 (("symamd: row %d col %d unsorted/duplicate\n", i, j)) ; |
|
|
1133 } |
|
|
1134 |
|
|
1135 if (i > j && mark [i] != j) |
|
|
1136 { |
|
|
1137 /* row k of M will contain column indices i and j */ |
|
|
1138 count [i]++ ; |
|
|
1139 count [j]++ ; |
|
|
1140 } |
|
|
1141 |
|
|
1142 /* mark the row as having been seen in this column */ |
|
|
1143 mark [i] = j ; |
|
|
1144 |
|
|
1145 last_row = i ; |
|
|
1146 } |
|
|
1147 } |
|
|
1148 |
|
|
1149 if (stats [COLAMD_STATUS] == COLAMD_OK) |
|
|
1150 { |
|
|
1151 /* if there are no duplicate entries, then mark is no longer needed */ |
|
|
1152 (*release) ((void *) mark) ; |
|
|
1153 } |
|
|
1154 |
|
|
1155 /* === Compute column pointers of M ===================================== */ |
|
|
1156 |
|
|
1157 /* use output permutation, perm, for column pointers of M */ |
|
|
1158 perm [0] = 0 ; |
|
|
1159 for (j = 1 ; j <= n ; j++) |
|
|
1160 { |
|
|
1161 perm [j] = perm [j-1] + count [j-1] ; |
|
|
1162 } |
|
|
1163 for (j = 0 ; j < n ; j++) |
|
|
1164 { |
|
|
1165 count [j] = perm [j] ; |
|
|
1166 } |
|
|
1167 |
|
|
1168 /* === Construct M ====================================================== */ |
|
|
1169 |
|
|
1170 mnz = perm [n] ; |
|
|
1171 n_row = mnz / 2 ; |
|
|
1172 Mlen = colamd_recommended (mnz, n_row, n) ; |
|
|
1173 M = (int *) ((*allocate) (Mlen, sizeof (int))) ; |
|
|
1174 DEBUG0 (("symamd: M is %d-by-%d with %d entries, Mlen = %d\n", |
|
|
1175 n_row, n, mnz, Mlen)) ; |
|
|
1176 |
|
|
1177 if (!M) |
|
|
1178 { |
|
|
1179 stats [COLAMD_STATUS] = COLAMD_ERROR_out_of_memory ; |
|
|
1180 (*release) ((void *) count) ; |
|
|
1181 (*release) ((void *) mark) ; |
|
|
1182 DEBUG0 (("symamd: allocate M (size %d) failed\n", Mlen)) ; |
|
|
1183 return (FALSE) ; |
|
|
1184 } |
|
|
1185 |
|
|
1186 k = 0 ; |
|
|
1187 |
|
|
1188 if (stats [COLAMD_STATUS] == COLAMD_OK) |
|
|
1189 { |
|
|
1190 /* Matrix is OK */ |
|
|
1191 for (j = 0 ; j < n ; j++) |
|
|
1192 { |
|
|
1193 ASSERT (p [j+1] - p [j] >= 0) ; |
|
|
1194 for (pp = p [j] ; pp < p [j+1] ; pp++) |
|
|
1195 { |
|
|
1196 i = A [pp] ; |
|
|
1197 ASSERT (i >= 0 && i < n) ; |
|
|
1198 if (i > j) |
|
|
1199 { |
|
|
1200 /* row k of M contains column indices i and j */ |
|
|
1201 M [count [i]++] = k ; |
|
|
1202 M [count [j]++] = k ; |
|
|
1203 k++ ; |
|
|
1204 } |
|
|
1205 } |
|
|
1206 } |
|
|
1207 } |
|
|
1208 else |
|
|
1209 { |
|
|
1210 /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */ |
|
|
1211 DEBUG0 (("symamd: Duplicates in A.\n")) ; |
|
|
1212 for (i = 0 ; i < n ; i++) |
|
|
1213 { |
|
|
1214 mark [i] = -1 ; |
|
|
1215 } |
|
|
1216 for (j = 0 ; j < n ; j++) |
|
|
1217 { |
|
|
1218 ASSERT (p [j+1] - p [j] >= 0) ; |
|
|
1219 for (pp = p [j] ; pp < p [j+1] ; pp++) |
|
|
1220 { |
|
|
1221 i = A [pp] ; |
|
|
1222 ASSERT (i >= 0 && i < n) ; |
|
|
1223 if (i > j && mark [i] != j) |
|
|
1224 { |
|
|
1225 /* row k of M contains column indices i and j */ |
|
|
1226 M [count [i]++] = k ; |
|
|
1227 M [count [j]++] = k ; |
|
|
1228 k++ ; |
|
|
1229 mark [i] = j ; |
|
|
1230 } |
|
|
1231 } |
|
|
1232 } |
|
|
1233 (*release) ((void *) mark) ; |
|
|
1234 } |
|
|
1235 |
|
|
1236 /* count and mark no longer needed */ |
|
|
1237 (*release) ((void *) count) ; |
|
|
1238 ASSERT (k == n_row) ; |
|
|
1239 |
|
|
1240 /* === Adjust the knobs for M =========================================== */ |
|
|
1241 |
|
|
1242 for (i = 0 ; i < COLAMD_KNOBS ; i++) |
|
|
1243 { |
|
|
1244 cknobs [i] = knobs [i] ; |
|
|
1245 } |
|
|
1246 |
|
|
1247 /* there are no dense rows in M */ |
|
|
1248 cknobs [COLAMD_DENSE_ROW] = 1.0 ; |
|
|
1249 |
|
|
1250 if (n_row != 0 && n < n_row) |
|
|
1251 { |
|
|
1252 /* On input, the knob is a fraction of 1..n, the number of rows of A. */ |
|
|
1253 /* Convert it to a fraction of 1..n_row, of the number of rows of M. */ |
|
|
1254 cknobs [COLAMD_DENSE_COL] = (knobs [COLAMD_DENSE_ROW] * n) / n_row ; |
|
|
1255 } |
|
|
1256 else |
|
|
1257 { |
|
|
1258 /* no dense columns in M */ |
|
|
1259 cknobs [COLAMD_DENSE_COL] = 1.0 ; |
|
|
1260 } |
|
|
1261 |
|
|
1262 DEBUG0 (("symamd: dense col knob for M: %g\n", cknobs [COLAMD_DENSE_COL])) ; |
|
|
1263 |
|
|
1264 /* === Order the columns of M =========================================== */ |
|
|
1265 |
|
|
1266 if (!colamd (n_row, n, Mlen, M, perm, cknobs, cstats)) |
|
|
1267 { |
|
|
1268 /* This "cannot" happen, unless there is a bug in the code. */ |
|
|
1269 stats [COLAMD_STATUS] = COLAMD_ERROR_internal_error ; |
|
|
1270 (*release) ((void *) M) ; |
|
|
1271 DEBUG0 (("symamd: internal error!\n")) ; |
|
|
1272 return (FALSE) ; |
|
|
1273 } |
|
|
1274 |
|
|
1275 /* Note that the output permutation is now in perm */ |
|
|
1276 |
|
|
1277 /* === get the statistics for symamd from colamd ======================== */ |
|
|
1278 |
|
|
1279 /* note that a dense column in colamd means a dense row and col in symamd */ |
|
|
1280 stats [COLAMD_DENSE_ROW] = cstats [COLAMD_DENSE_COL] ; |
|
|
1281 stats [COLAMD_DENSE_COL] = cstats [COLAMD_DENSE_COL] ; |
|
|
1282 stats [COLAMD_DEFRAG_COUNT] = cstats [COLAMD_DEFRAG_COUNT] ; |
|
|
1283 |
|
|
1284 /* === Free M =========================================================== */ |
|
|
1285 |
|
|
1286 (*release) ((void *) M) ; |
|
|
1287 DEBUG0 (("symamd: done.\n")) ; |
|
|
1288 return (TRUE) ; |
|
|
1289 |
|
|
1290 } |
|
|
1291 |
|
|
1292 /* ========================================================================== */ |
|
|
1293 /* === colamd =============================================================== */ |
|
|
1294 /* ========================================================================== */ |
|
|
1295 |
|
|
1296 /* |
|
|
1297 The colamd routine computes a column ordering Q of a sparse matrix |
|
|
1298 A such that the LU factorization P(AQ) = LU remains sparse, where P is |
|
|
1299 selected via partial pivoting. The routine can also be viewed as |
|
|
1300 providing a permutation Q such that the Cholesky factorization |
|
|
1301 (AQ)'(AQ) = LL' remains sparse. |
|
|
1302 */ |
|
|
1303 |
|
|
1304 PUBLIC int colamd /* returns TRUE if successful, FALSE otherwise*/ |
|
|
1305 ( |
|
|
1306 /* === Parameters ======================================================= */ |
|
|
1307 |
|
|
1308 int n_row, /* number of rows in A */ |
|
|
1309 int n_col, /* number of columns in A */ |
|
|
1310 int Alen, /* length of A */ |
|
|
1311 int A [], /* row indices of A */ |
|
|
1312 int p [], /* pointers to columns in A */ |
|
|
1313 double knobs [COLAMD_KNOBS],/* parameters (uses defaults if NULL) */ |
|
|
1314 int stats [COLAMD_STATS] /* output statistics and error codes */ |
|
|
1315 ) |
|
|
1316 { |
|
|
1317 /* === Local variables ================================================== */ |
|
|
1318 |
|
|
1319 int i ; /* loop index */ |
|
|
1320 int nnz ; /* nonzeros in A */ |
|
|
1321 int Row_size ; /* size of Row [], in integers */ |
|
|
1322 int Col_size ; /* size of Col [], in integers */ |
|
|
1323 int need ; /* minimum required length of A */ |
|
|
1324 Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ |
|
|
1325 Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ |
|
|
1326 int n_col2 ; /* number of non-dense, non-empty columns */ |
|
|
1327 int n_row2 ; /* number of non-dense, non-empty rows */ |
|
|
1328 int ngarbage ; /* number of garbage collections performed */ |
|
|
1329 int max_deg ; /* maximum row degree */ |
|
|
1330 double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ |
|
|
1331 |
|
|
1332 #ifndef NDEBUG |
|
|
1333 colamd_get_debug ("colamd") ; |
|
|
1334 #endif /* NDEBUG */ |
|
|
1335 |
|
|
1336 /* === Check the input arguments ======================================== */ |
|
|
1337 |
|
|
1338 if (!stats) |
|
|
1339 { |
|
|
1340 DEBUG0 (("colamd: stats not present\n")) ; |
|
|
1341 return (FALSE) ; |
|
|
1342 } |
|
|
1343 for (i = 0 ; i < COLAMD_STATS ; i++) |
|
|
1344 { |
|
|
1345 stats [i] = 0 ; |
|
|
1346 } |
|
|
1347 stats [COLAMD_STATUS] = COLAMD_OK ; |
|
|
1348 stats [COLAMD_INFO1] = -1 ; |
|
|
1349 stats [COLAMD_INFO2] = -1 ; |
|
|
1350 |
|
|
1351 if (!A) /* A is not present */ |
|
|
1352 { |
|
|
1353 stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; |
|
|
1354 DEBUG0 (("colamd: A not present\n")) ; |
|
|
1355 return (FALSE) ; |
|
|
1356 } |
|
|
1357 |
|
|
1358 if (!p) /* p is not present */ |
|
|
1359 { |
|
|
1360 stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; |
|
|
1361 DEBUG0 (("colamd: p not present\n")) ; |
|
|
1362 return (FALSE) ; |
|
|
1363 } |
|
|
1364 |
|
|
1365 if (n_row < 0) /* n_row must be >= 0 */ |
|
|
1366 { |
|
|
1367 stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; |
|
|
1368 stats [COLAMD_INFO1] = n_row ; |
|
|
1369 DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; |
|
|
1370 return (FALSE) ; |
|
|
1371 } |
|
|
1372 |
|
|
1373 if (n_col < 0) /* n_col must be >= 0 */ |
|
|
1374 { |
|
|
1375 stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; |
|
|
1376 stats [COLAMD_INFO1] = n_col ; |
|
|
1377 DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; |
|
|
1378 return (FALSE) ; |
|
|
1379 } |
|
|
1380 |
|
|
1381 nnz = p [n_col] ; |
|
|
1382 if (nnz < 0) /* nnz must be >= 0 */ |
|
|
1383 { |
|
|
1384 stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; |
|
|
1385 stats [COLAMD_INFO1] = nnz ; |
|
|
1386 DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; |
|
|
1387 return (FALSE) ; |
|
|
1388 } |
|
|
1389 |
|
|
1390 if (p [0] != 0) |
|
|
1391 { |
|
|
1392 stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; |
|
|
1393 stats [COLAMD_INFO1] = p [0] ; |
|
|
1394 DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; |
|
|
1395 return (FALSE) ; |
|
|
1396 } |
|
|
1397 |
|
|
1398 /* === If no knobs, set default knobs =================================== */ |
|
|
1399 |
|
|
1400 if (!knobs) |
|
|
1401 { |
|
|
1402 colamd_set_defaults (default_knobs) ; |
|
|
1403 knobs = default_knobs ; |
|
|
1404 } |
|
|
1405 |
|
|
1406 /* === Allocate the Row and Col arrays from array A ===================== */ |
|
|
1407 |
|
|
1408 Col_size = COLAMD_C (n_col) ; |
|
|
1409 Row_size = COLAMD_R (n_row) ; |
|
|
1410 need = 2*nnz + n_col + Col_size + Row_size ; |
|
|
1411 |
|
|
1412 if (need > Alen) |
|
|
1413 { |
|
|
1414 /* not enough space in array A to perform the ordering */ |
|
|
1415 stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; |
|
|
1416 stats [COLAMD_INFO1] = need ; |
|
|
1417 stats [COLAMD_INFO2] = Alen ; |
|
|
1418 DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); |
|
|
1419 return (FALSE) ; |
|
|
1420 } |
|
|
1421 |
|
|
1422 Alen -= Col_size + Row_size ; |
|
|
1423 Col = (Colamd_Col *) &A [Alen] ; |
|
|
1424 Row = (Colamd_Row *) &A [Alen + Col_size] ; |
|
|
1425 |
|
|
1426 /* === Construct the row and column data structures ===================== */ |
|
|
1427 |
|
|
1428 if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) |
|
|
1429 { |
|
|
1430 /* input matrix is invalid */ |
|
|
1431 DEBUG0 (("colamd: Matrix invalid\n")) ; |
|
|
1432 return (FALSE) ; |
|
|
1433 } |
|
|
1434 |
|
|
1435 /* === Initialize scores, kill dense rows/columns ======================= */ |
|
|
1436 |
|
|
1437 init_scoring (n_row, n_col, Row, Col, A, p, knobs, |
|
|
1438 &n_row2, &n_col2, &max_deg) ; |
|
|
1439 |
|
|
1440 /* === Order the supercolumns =========================================== */ |
|
|
1441 |
|
|
1442 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, |
|
|
1443 n_col2, max_deg, 2*nnz) ; |
|
|
1444 |
|
|
1445 /* === Order the non-principal columns ================================== */ |
|
|
1446 |
|
|
1447 order_children (n_col, Col, p) ; |
|
|
1448 |
|
|
1449 /* === Return statistics in stats ======================================= */ |
|
|
1450 |
|
|
1451 stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; |
|
|
1452 stats [COLAMD_DENSE_COL] = n_col - n_col2 ; |
|
|
1453 stats [COLAMD_DEFRAG_COUNT] = ngarbage ; |
|
|
1454 DEBUG0 (("colamd: done.\n")) ; |
|
|
1455 return (TRUE) ; |
|
|
1456 } |
|
|
1457 |
|
|
1458 |
|
|
1459 /* ========================================================================== */ |
|
|
1460 /* === colamd_report ======================================================== */ |
|
|
1461 /* ========================================================================== */ |
|
|
1462 |
|
|
1463 PUBLIC void colamd_report |
|
|
1464 ( |
|
|
1465 int stats [COLAMD_STATS] |
|
|
1466 ) |
|
|
1467 { |
|
|
1468 print_report ("colamd", stats) ; |
|
|
1469 } |
|
|
1470 |
|
|
1471 |
|
|
1472 /* ========================================================================== */ |
|
|
1473 /* === symamd_report ======================================================== */ |
|
|
1474 /* ========================================================================== */ |
|
|
1475 |
|
|
1476 PUBLIC void symamd_report |
|
|
1477 ( |
|
|
1478 int stats [COLAMD_STATS] |
|
|
1479 ) |
|
|
1480 { |
|
|
1481 print_report ("symamd", stats) ; |
|
|
1482 } |
|
|
1483 |
|
|
1484 |
|
|
1485 |
|
|
1486 /* ========================================================================== */ |
|
|
1487 /* === NON-USER-CALLABLE ROUTINES: ========================================== */ |
|
|
1488 /* ========================================================================== */ |
|
|
1489 |
|
|
1490 /* There are no user-callable routines beyond this point in the file */ |
|
|
1491 |
|
|
1492 |
|
|
1493 /* ========================================================================== */ |
|
|
1494 /* === init_rows_cols ======================================================= */ |
|
|
1495 /* ========================================================================== */ |
|
|
1496 |
|
|
1497 /* |
|
|
1498 Takes the column form of the matrix in A and creates the row form of the |
|
|
1499 matrix. Also, row and column attributes are stored in the Col and Row |
|
|
1500 structs. If the columns are un-sorted or contain duplicate row indices, |
|
|
1501 this routine will also sort and remove duplicate row indices from the |
|
|
1502 column form of the matrix. Returns FALSE if the matrix is invalid, |
|
|
1503 TRUE otherwise. Not user-callable. |
|
|
1504 */ |
|
|
1505 |
|
|
1506 PRIVATE int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */ |
|
|
1507 ( |
|
|
1508 /* === Parameters ======================================================= */ |
|
|
1509 |
|
|
1510 int n_row, /* number of rows of A */ |
|
|
1511 int n_col, /* number of columns of A */ |
|
|
1512 Colamd_Row Row [], /* of size n_row+1 */ |
|
|
1513 Colamd_Col Col [], /* of size n_col+1 */ |
|
|
1514 int A [], /* row indices of A, of size Alen */ |
|
|
1515 int p [], /* pointers to columns in A, of size n_col+1 */ |
|
|
1516 int stats [COLAMD_STATS] /* colamd statistics */ |
|
|
1517 ) |
|
|
1518 { |
|
|
1519 /* === Local variables ================================================== */ |
|
|
1520 |
|
|
1521 int col ; /* a column index */ |
|
|
1522 int row ; /* a row index */ |
|
|
1523 int *cp ; /* a column pointer */ |
|
|
1524 int *cp_end ; /* a pointer to the end of a column */ |
|
|
1525 int *rp ; /* a row pointer */ |
|
|
1526 int *rp_end ; /* a pointer to the end of a row */ |
|
|
1527 int last_row ; /* previous row */ |
|
|
1528 |
|
|
1529 /* === Initialize columns, and check column pointers ==================== */ |
|
|
1530 |
|
|
1531 for (col = 0 ; col < n_col ; col++) |
|
|
1532 { |
|
|
1533 Col [col].start = p [col] ; |
|
|
1534 Col [col].length = p [col+1] - p [col] ; |
|
|
1535 |
|
|
1536 if (Col [col].length < 0) |
|
|
1537 { |
|
|
1538 /* column pointers must be non-decreasing */ |
|
|
1539 stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; |
|
|
1540 stats [COLAMD_INFO1] = col ; |
|
|
1541 stats [COLAMD_INFO2] = Col [col].length ; |
|
|
1542 DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; |
|
|
1543 return (FALSE) ; |
|
|
1544 } |
|
|
1545 |
|
|
1546 Col [col].shared1.thickness = 1 ; |
|
|
1547 Col [col].shared2.score = 0 ; |
|
|
1548 Col [col].shared3.prev = EMPTY ; |
|
|
1549 Col [col].shared4.degree_next = EMPTY ; |
|
|
1550 } |
|
|
1551 |
|
|
1552 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ |
|
|
1553 |
|
|
1554 /* === Scan columns, compute row degrees, and check row indices ========= */ |
|
|
1555 |
|
|
1556 stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ |
|
|
1557 |
|
|
1558 for (row = 0 ; row < n_row ; row++) |
|
|
1559 { |
|
|
1560 Row [row].length = 0 ; |
|
|
1561 Row [row].shared2.mark = -1 ; |
|
|
1562 } |
|
|
1563 |
|
|
1564 for (col = 0 ; col < n_col ; col++) |
|
|
1565 { |
|
|
1566 last_row = -1 ; |
|
|
1567 |
|
|
1568 cp = &A [p [col]] ; |
|
|
1569 cp_end = &A [p [col+1]] ; |
|
|
1570 |
|
|
1571 while (cp < cp_end) |
|
|
1572 { |
|
|
1573 row = *cp++ ; |
|
|
1574 |
|
|
1575 /* make sure row indices within range */ |
|
|
1576 if (row < 0 || row >= n_row) |
|
|
1577 { |
|
|
1578 stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; |
|
|
1579 stats [COLAMD_INFO1] = col ; |
|
|
1580 stats [COLAMD_INFO2] = row ; |
|
|
1581 stats [COLAMD_INFO3] = n_row ; |
|
|
1582 DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; |
|
|
1583 return (FALSE) ; |
|
|
1584 } |
|
|
1585 |
|
|
1586 if (row <= last_row || Row [row].shared2.mark == col) |
|
|
1587 { |
|
|
1588 /* row index are unsorted or repeated (or both), thus col */ |
|
|
1589 /* is jumbled. This is a notice, not an error condition. */ |
|
|
1590 stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; |
|
|
1591 stats [COLAMD_INFO1] = col ; |
|
|
1592 stats [COLAMD_INFO2] = row ; |
|
|
1593 (stats [COLAMD_INFO3]) ++ ; |
|
|
1594 DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); |
|
|
1595 } |
|
|
1596 |
|
|
1597 if (Row [row].shared2.mark != col) |
|
|
1598 { |
|
|
1599 Row [row].length++ ; |
|
|
1600 } |
|
|
1601 else |
|
|
1602 { |
|
|
1603 /* this is a repeated entry in the column, */ |
|
|
1604 /* it will be removed */ |
|
|
1605 Col [col].length-- ; |
|
|
1606 } |
|
|
1607 |
|
|
1608 /* mark the row as having been seen in this column */ |
|
|
1609 Row [row].shared2.mark = col ; |
|
|
1610 |
|
|
1611 last_row = row ; |
|
|
1612 } |
|
|
1613 } |
|
|
1614 |
|
|
1615 /* === Compute row pointers ============================================= */ |
|
|
1616 |
|
|
1617 /* row form of the matrix starts directly after the column */ |
|
|
1618 /* form of matrix in A */ |
|
|
1619 Row [0].start = p [n_col] ; |
|
|
1620 Row [0].shared1.p = Row [0].start ; |
|
|
1621 Row [0].shared2.mark = -1 ; |
|
|
1622 for (row = 1 ; row < n_row ; row++) |
|
|
1623 { |
|
|
1624 Row [row].start = Row [row-1].start + Row [row-1].length ; |
|
|
1625 Row [row].shared1.p = Row [row].start ; |
|
|
1626 Row [row].shared2.mark = -1 ; |
|
|
1627 } |
|
|
1628 |
|
|
1629 /* === Create row form ================================================== */ |
|
|
1630 |
|
|
1631 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) |
|
|
1632 { |
|
|
1633 /* if cols jumbled, watch for repeated row indices */ |
|
|
1634 for (col = 0 ; col < n_col ; col++) |
|
|
1635 { |
|
|
1636 cp = &A [p [col]] ; |
|
|
1637 cp_end = &A [p [col+1]] ; |
|
|
1638 while (cp < cp_end) |
|
|
1639 { |
|
|
1640 row = *cp++ ; |
|
|
1641 if (Row [row].shared2.mark != col) |
|
|
1642 { |
|
|
1643 A [(Row [row].shared1.p)++] = col ; |
|
|
1644 Row [row].shared2.mark = col ; |
|
|
1645 } |
|
|
1646 } |
|
|
1647 } |
|
|
1648 } |
|
|
1649 else |
|
|
1650 { |
|
|
1651 /* if cols not jumbled, we don't need the mark (this is faster) */ |
|
|
1652 for (col = 0 ; col < n_col ; col++) |
|
|
1653 { |
|
|
1654 cp = &A [p [col]] ; |
|
|
1655 cp_end = &A [p [col+1]] ; |
|
|
1656 while (cp < cp_end) |
|
|
1657 { |
|
|
1658 A [(Row [*cp++].shared1.p)++] = col ; |
|
|
1659 } |
|
|
1660 } |
|
|
1661 } |
|
|
1662 |
|
|
1663 /* === Clear the row marks and set row degrees ========================== */ |
|
|
1664 |
|
|
1665 for (row = 0 ; row < n_row ; row++) |
|
|
1666 { |
|
|
1667 Row [row].shared2.mark = 0 ; |
|
|
1668 Row [row].shared1.degree = Row [row].length ; |
|
|
1669 } |
|
|
1670 |
|
|
1671 /* === See if we need to re-create columns ============================== */ |
|
|
1672 |
|
|
1673 if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) |
|
|
1674 { |
|
|
1675 DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; |
|
|
1676 |
|
|
1677 #ifndef NDEBUG |
|
|
1678 /* make sure column lengths are correct */ |
|
|
1679 for (col = 0 ; col < n_col ; col++) |
|
|
1680 { |
|
|
1681 p [col] = Col [col].length ; |
|
|
1682 } |
|
|
1683 for (row = 0 ; row < n_row ; row++) |
|
|
1684 { |
|
|
1685 rp = &A [Row [row].start] ; |
|
|
1686 rp_end = rp + Row [row].length ; |
|
|
1687 while (rp < rp_end) |
|
|
1688 { |
|
|
1689 p [*rp++]-- ; |
|
|
1690 } |
|
|
1691 } |
|
|
1692 for (col = 0 ; col < n_col ; col++) |
|
|
1693 { |
|
|
1694 ASSERT (p [col] == 0) ; |
|
|
1695 } |
|
|
1696 /* now p is all zero (different than when debugging is turned off) */ |
|
|
1697 #endif /* NDEBUG */ |
|
|
1698 |
|
|
1699 /* === Compute col pointers ========================================= */ |
|
|
1700 |
|
|
1701 /* col form of the matrix starts at A [0]. */ |
|
|
1702 /* Note, we may have a gap between the col form and the row */ |
|
|
1703 /* form if there were duplicate entries, if so, it will be */ |
|
|
1704 /* removed upon the first garbage collection */ |
|
|
1705 Col [0].start = 0 ; |
|
|
1706 p [0] = Col [0].start ; |
|
|
1707 for (col = 1 ; col < n_col ; col++) |
|
|
1708 { |
|
|
1709 /* note that the lengths here are for pruned columns, i.e. */ |
|
|
1710 /* no duplicate row indices will exist for these columns */ |
|
|
1711 Col [col].start = Col [col-1].start + Col [col-1].length ; |
|
|
1712 p [col] = Col [col].start ; |
|
|
1713 } |
|
|
1714 |
|
|
1715 /* === Re-create col form =========================================== */ |
|
|
1716 |
|
|
1717 for (row = 0 ; row < n_row ; row++) |
|
|
1718 { |
|
|
1719 rp = &A [Row [row].start] ; |
|
|
1720 rp_end = rp + Row [row].length ; |
|
|
1721 while (rp < rp_end) |
|
|
1722 { |
|
|
1723 A [(p [*rp++])++] = row ; |
|
|
1724 } |
|
|
1725 } |
|
|
1726 } |
|
|
1727 |
|
|
1728 /* === Done. Matrix is not (or no longer) jumbled ====================== */ |
|
|
1729 |
|
|
1730 return (TRUE) ; |
|
|
1731 } |
|
|
1732 |
|
|
1733 |
|
|
1734 /* ========================================================================== */ |
|
|
1735 /* === init_scoring ========================================================= */ |
|
|
1736 /* ========================================================================== */ |
|
|
1737 |
|
|
1738 /* |
|
|
1739 Kills dense or empty columns and rows, calculates an initial score for |
|
|
1740 each column, and places all columns in the degree lists. Not user-callable. |
|
|
1741 */ |
|
|
1742 |
|
|
1743 PRIVATE void init_scoring |
|
|
1744 ( |
|
|
1745 /* === Parameters ======================================================= */ |
|
|
1746 |
|
|
1747 int n_row, /* number of rows of A */ |
|
|
1748 int n_col, /* number of columns of A */ |
|
|
1749 Colamd_Row Row [], /* of size n_row+1 */ |
|
|
1750 Colamd_Col Col [], /* of size n_col+1 */ |
|
|
1751 int A [], /* column form and row form of A */ |
|
|
1752 int head [], /* of size n_col+1 */ |
|
|
1753 double knobs [COLAMD_KNOBS],/* parameters */ |
|
|
1754 int *p_n_row2, /* number of non-dense, non-empty rows */ |
|
|
1755 int *p_n_col2, /* number of non-dense, non-empty columns */ |
|
|
1756 int *p_max_deg /* maximum row degree */ |
|
|
1757 ) |
|
|
1758 { |
|
|
1759 /* === Local variables ================================================== */ |
|
|
1760 |
|
|
1761 int c ; /* a column index */ |
|
|
1762 int r, row ; /* a row index */ |
|
|
1763 int *cp ; /* a column pointer */ |
|
|
1764 int deg ; /* degree of a row or column */ |
|
|
1765 int *cp_end ; /* a pointer to the end of a column */ |
|
|
1766 int *new_cp ; /* new column pointer */ |
|
|
1767 int col_length ; /* length of pruned column */ |
|
|
1768 int score ; /* current column score */ |
|
|
1769 int n_col2 ; /* number of non-dense, non-empty columns */ |
|
|
1770 int n_row2 ; /* number of non-dense, non-empty rows */ |
|
|
1771 int dense_row_count ; /* remove rows with more entries than this */ |
|
|
1772 int dense_col_count ; /* remove cols with more entries than this */ |
|
|
1773 int min_score ; /* smallest column score */ |
|
|
1774 int max_deg ; /* maximum row degree */ |
|
|
1775 int next_col ; /* Used to add to degree list.*/ |
|
|
1776 |
|
|
1777 #ifndef NDEBUG |
|
|
1778 int debug_count ; /* debug only. */ |
|
|
1779 #endif /* NDEBUG */ |
|
|
1780 |
|
|
1781 /* === Extract knobs ==================================================== */ |
|
|
1782 |
|
|
1783 dense_row_count = MAX (0, MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ; |
|
|
1784 dense_col_count = MAX (0, MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ; |
|
|
1785 DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; |
|
|
1786 max_deg = 0 ; |
|
|
1787 n_col2 = n_col ; |
|
|
1788 n_row2 = n_row ; |
|
|
1789 |
|
|
1790 /* === Kill empty columns =============================================== */ |
|
|
1791 |
|
|
1792 /* Put the empty columns at the end in their natural order, so that LU */ |
|
|
1793 /* factorization can proceed as far as possible. */ |
|
|
1794 for (c = n_col-1 ; c >= 0 ; c--) |
|
|
1795 { |
|
|
1796 deg = Col [c].length ; |
|
|
1797 if (deg == 0) |
|
|
1798 { |
|
|
1799 /* this is a empty column, kill and order it last */ |
|
|
1800 Col [c].shared2.order = --n_col2 ; |
|
|
1801 KILL_PRINCIPAL_COL (c) ; |
|
|
1802 } |
|
|
1803 } |
|
|
1804 DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; |
|
|
1805 |
|
|
1806 /* === Kill dense columns =============================================== */ |
|
|
1807 |
|
|
1808 /* Put the dense columns at the end, in their natural order */ |
|
|
1809 for (c = n_col-1 ; c >= 0 ; c--) |
|
|
1810 { |
|
|
1811 /* skip any dead columns */ |
|
|
1812 if (COL_IS_DEAD (c)) |
|
|
1813 { |
|
|
1814 continue ; |
|
|
1815 } |
|
|
1816 deg = Col [c].length ; |
|
|
1817 if (deg > dense_col_count) |
|
|
1818 { |
|
|
1819 /* this is a dense column, kill and order it last */ |
|
|
1820 Col [c].shared2.order = --n_col2 ; |
|
|
1821 /* decrement the row degrees */ |
|
|
1822 cp = &A [Col [c].start] ; |
|
|
1823 cp_end = cp + Col [c].length ; |
|
|
1824 while (cp < cp_end) |
|
|
1825 { |
|
|
1826 Row [*cp++].shared1.degree-- ; |
|
|
1827 } |
|
|
1828 KILL_PRINCIPAL_COL (c) ; |
|
|
1829 } |
|
|
1830 } |
|
|
1831 DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; |
|
|
1832 |
|
|
1833 /* === Kill dense and empty rows ======================================== */ |
|
|
1834 |
|
|
1835 for (r = 0 ; r < n_row ; r++) |
|
|
1836 { |
|
|
1837 deg = Row [r].shared1.degree ; |
|
|
1838 ASSERT (deg >= 0 && deg <= n_col) ; |
|
|
1839 if (deg > dense_row_count || deg == 0) |
|
|
1840 { |
|
|
1841 /* kill a dense or empty row */ |
|
|
1842 KILL_ROW (r) ; |
|
|
1843 --n_row2 ; |
|
|
1844 } |
|
|
1845 else |
|
|
1846 { |
|
|
1847 /* keep track of max degree of remaining rows */ |
|
|
1848 max_deg = MAX (max_deg, deg) ; |
|
|
1849 } |
|
|
1850 } |
|
|
1851 DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; |
|
|
1852 |
|
|
1853 /* === Compute initial column scores ==================================== */ |
|
|
1854 |
|
|
1855 /* At this point the row degrees are accurate. They reflect the number */ |
|
|
1856 /* of "live" (non-dense) columns in each row. No empty rows exist. */ |
|
|
1857 /* Some "live" columns may contain only dead rows, however. These are */ |
|
|
1858 /* pruned in the code below. */ |
|
|
1859 |
|
|
1860 /* now find the initial matlab score for each column */ |
|
|
1861 for (c = n_col-1 ; c >= 0 ; c--) |
|
|
1862 { |
|
|
1863 /* skip dead column */ |
|
|
1864 if (COL_IS_DEAD (c)) |
|
|
1865 { |
|
|
1866 continue ; |
|
|
1867 } |
|
|
1868 score = 0 ; |
|
|
1869 cp = &A [Col [c].start] ; |
|
|
1870 new_cp = cp ; |
|
|
1871 cp_end = cp + Col [c].length ; |
|
|
1872 while (cp < cp_end) |
|
|
1873 { |
|
|
1874 /* get a row */ |
|
|
1875 row = *cp++ ; |
|
|
1876 /* skip if dead */ |
|
|
1877 if (ROW_IS_DEAD (row)) |
|
|
1878 { |
|
|
1879 continue ; |
|
|
1880 } |
|
|
1881 /* compact the column */ |
|
|
1882 *new_cp++ = row ; |
|
|
1883 /* add row's external degree */ |
|
|
1884 score += Row [row].shared1.degree - 1 ; |
|
|
1885 /* guard against integer overflow */ |
|
|
1886 score = MIN (score, n_col) ; |
|
|
1887 } |
|
|
1888 /* determine pruned column length */ |
|
|
1889 col_length = (int) (new_cp - &A [Col [c].start]) ; |
|
|
1890 if (col_length == 0) |
|
|
1891 { |
|
|
1892 /* a newly-made null column (all rows in this col are "dense" */ |
|
|
1893 /* and have already been killed) */ |
|
|
1894 DEBUG2 (("Newly null killed: %d\n", c)) ; |
|
|
1895 Col [c].shared2.order = --n_col2 ; |
|
|
1896 KILL_PRINCIPAL_COL (c) ; |
|
|
1897 } |
|
|
1898 else |
|
|
1899 { |
|
|
1900 /* set column length and set score */ |
|
|
1901 ASSERT (score >= 0) ; |
|
|
1902 ASSERT (score <= n_col) ; |
|
|
1903 Col [c].length = col_length ; |
|
|
1904 Col [c].shared2.score = score ; |
|
|
1905 } |
|
|
1906 } |
|
|
1907 DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", |
|
|
1908 n_col-n_col2)) ; |
|
|
1909 |
|
|
1910 /* At this point, all empty rows and columns are dead. All live columns */ |
|
|
1911 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */ |
|
|
1912 /* yet). Rows may contain dead columns, but all live rows contain at */ |
|
|
1913 /* least one live column. */ |
|
|
1914 |
|
|
1915 #ifndef NDEBUG |
|
|
1916 debug_structures (n_row, n_col, Row, Col, A, n_col2) ; |
|
|
1917 #endif /* NDEBUG */ |
|
|
1918 |
|
|
1919 /* === Initialize degree lists ========================================== */ |
|
|
1920 |
|
|
1921 #ifndef NDEBUG |
|
|
1922 debug_count = 0 ; |
|
|
1923 #endif /* NDEBUG */ |
|
|
1924 |
|
|
1925 /* clear the hash buckets */ |
|
|
1926 for (c = 0 ; c <= n_col ; c++) |
|
|
1927 { |
|
|
1928 head [c] = EMPTY ; |
|
|
1929 } |
|
|
1930 min_score = n_col ; |
|
|
1931 /* place in reverse order, so low column indices are at the front */ |
|
|
1932 /* of the lists. This is to encourage natural tie-breaking */ |
|
|
1933 for (c = n_col-1 ; c >= 0 ; c--) |
|
|
1934 { |
|
|
1935 /* only add principal columns to degree lists */ |
|
|
1936 if (COL_IS_ALIVE (c)) |
|
|
1937 { |
|
|
1938 DEBUG4 (("place %d score %d minscore %d ncol %d\n", |
|
|
1939 c, Col [c].shared2.score, min_score, n_col)) ; |
|
|
1940 |
|
|
1941 /* === Add columns score to DList =============================== */ |
|
|
1942 |
|
|
1943 score = Col [c].shared2.score ; |
|
|
1944 |
|
|
1945 ASSERT (min_score >= 0) ; |
|
|
1946 ASSERT (min_score <= n_col) ; |
|
|
1947 ASSERT (score >= 0) ; |
|
|
1948 ASSERT (score <= n_col) ; |
|
|
1949 ASSERT (head [score] >= EMPTY) ; |
|
|
1950 |
|
|
1951 /* now add this column to dList at proper score location */ |
|
|
1952 next_col = head [score] ; |
|
|
1953 Col [c].shared3.prev = EMPTY ; |
|
|
1954 Col [c].shared4.degree_next = next_col ; |
|
|
1955 |
|
|
1956 /* if there already was a column with the same score, set its */ |
|
|
1957 /* previous pointer to this new column */ |
|
|
1958 if (next_col != EMPTY) |
|
|
1959 { |
|
|
1960 Col [next_col].shared3.prev = c ; |
|
|
1961 } |
|
|
1962 head [score] = c ; |
|
|
1963 |
|
|
1964 /* see if this score is less than current min */ |
|
|
1965 min_score = MIN (min_score, score) ; |
|
|
1966 |
|
|
1967 #ifndef NDEBUG |
|
|
1968 debug_count++ ; |
|
|
1969 #endif /* NDEBUG */ |
|
|
1970 |
|
|
1971 } |
|
|
1972 } |
|
|
1973 |
|
|
1974 #ifndef NDEBUG |
|
|
1975 DEBUG1 (("colamd: Live cols %d out of %d, non-princ: %d\n", |
|
|
1976 debug_count, n_col, n_col-debug_count)) ; |
|
|
1977 ASSERT (debug_count == n_col2) ; |
|
|
1978 debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; |
|
|
1979 #endif /* NDEBUG */ |
|
|
1980 |
|
|
1981 /* === Return number of remaining columns, and max row degree =========== */ |
|
|
1982 |
|
|
1983 *p_n_col2 = n_col2 ; |
|
|
1984 *p_n_row2 = n_row2 ; |
|
|
1985 *p_max_deg = max_deg ; |
|
|
1986 } |
|
|
1987 |
|
|
1988 |
|
|
1989 /* ========================================================================== */ |
|
|
1990 /* === find_ordering ======================================================== */ |
|
|
1991 /* ========================================================================== */ |
|
|
1992 |
|
|
1993 /* |
|
|
1994 Order the principal columns of the supercolumn form of the matrix |
|
|
1995 (no supercolumns on input). Uses a minimum approximate column minimum |
|
|
1996 degree ordering method. Not user-callable. |
|
|
1997 */ |
|
|
1998 |
|
|
1999 PRIVATE int find_ordering /* return the number of garbage collections */ |
|
|
2000 ( |
|
|
2001 /* === Parameters ======================================================= */ |
|
|
2002 |
|
|
2003 int n_row, /* number of rows of A */ |
|
|
2004 int n_col, /* number of columns of A */ |
|
|
2005 int Alen, /* size of A, 2*nnz + n_col or larger */ |
|
|
2006 Colamd_Row Row [], /* of size n_row+1 */ |
|
|
2007 Colamd_Col Col [], /* of size n_col+1 */ |
|
|
2008 int A [], /* column form and row form of A */ |
|
|
2009 int head [], /* of size n_col+1 */ |
|
|
2010 int n_col2, /* Remaining columns to order */ |
|
|
2011 int max_deg, /* Maximum row degree */ |
|
|
2012 int pfree /* index of first free slot (2*nnz on entry) */ |
|
|
2013 ) |
|
|
2014 { |
|
|
2015 /* === Local variables ================================================== */ |
|
|
2016 |
|
|
2017 int k ; /* current pivot ordering step */ |
|
|
2018 int pivot_col ; /* current pivot column */ |
|
|
2019 int *cp ; /* a column pointer */ |
|
|
2020 int *rp ; /* a row pointer */ |
|
|
2021 int pivot_row ; /* current pivot row */ |
|
|
2022 int *new_cp ; /* modified column pointer */ |
|
|
2023 int *new_rp ; /* modified row pointer */ |
|
|
2024 int pivot_row_start ; /* pointer to start of pivot row */ |
|
|
2025 int pivot_row_degree ; /* number of columns in pivot row */ |
|
|
2026 int pivot_row_length ; /* number of supercolumns in pivot row */ |
|
|
2027 int pivot_col_score ; /* score of pivot column */ |
|
|
2028 int needed_memory ; /* free space needed for pivot row */ |
|
|
2029 int *cp_end ; /* pointer to the end of a column */ |
|
|
2030 int *rp_end ; /* pointer to the end of a row */ |
|
|
2031 int row ; /* a row index */ |
|
|
2032 int col ; /* a column index */ |
|
|
2033 int max_score ; /* maximum possible score */ |
|
|
2034 int cur_score ; /* score of current column */ |
|
|
2035 unsigned int hash ; /* hash value for supernode detection */ |
|
|
2036 int head_column ; /* head of hash bucket */ |
|
|
2037 int first_col ; /* first column in hash bucket */ |
|
|
2038 int tag_mark ; /* marker value for mark array */ |
|
|
2039 int row_mark ; /* Row [row].shared2.mark */ |
|
|
2040 int set_difference ; /* set difference size of row with pivot row */ |
|
|
2041 int min_score ; /* smallest column score */ |
|
|
2042 int col_thickness ; /* "thickness" (no. of columns in a supercol) */ |
|
|
2043 int max_mark ; /* maximum value of tag_mark */ |
|
|
2044 int pivot_col_thickness ; /* number of columns represented by pivot col */ |
|
|
2045 int prev_col ; /* Used by Dlist operations. */ |
|
|
2046 int next_col ; /* Used by Dlist operations. */ |
|
|
2047 int ngarbage ; /* number of garbage collections performed */ |
|
|
2048 |
|
|
2049 #ifndef NDEBUG |
|
|
2050 int debug_d ; /* debug loop counter */ |
|
|
2051 int debug_step = 0 ; /* debug loop counter */ |
|
|
2052 #endif /* NDEBUG */ |
|
|
2053 |
|
|
2054 /* === Initialization and clear mark ==================================== */ |
|
|
2055 |
|
|
2056 max_mark = INT_MAX - n_col ; /* INT_MAX defined in <limits.h> */ |
|
|
2057 tag_mark = clear_mark (n_row, Row) ; |
|
|
2058 min_score = 0 ; |
|
|
2059 ngarbage = 0 ; |
|
|
2060 DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; |
|
|
2061 |
|
|
2062 /* === Order the columns ================================================ */ |
|
|
2063 |
|
|
2064 for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) |
|
|
2065 { |
|
|
2066 |
|
|
2067 #ifndef NDEBUG |
|
|
2068 if (debug_step % 100 == 0) |
|
|
2069 { |
|
|
2070 DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; |
|
|
2071 } |
|
|
2072 else |
|
|
2073 { |
|
|
2074 DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; |
|
|
2075 } |
|
|
2076 debug_step++ ; |
|
|
2077 debug_deg_lists (n_row, n_col, Row, Col, head, |
|
|
2078 min_score, n_col2-k, max_deg) ; |
|
|
2079 debug_matrix (n_row, n_col, Row, Col, A) ; |
|
|
2080 #endif /* NDEBUG */ |
|
|
2081 |
|
|
2082 /* === Select pivot column, and order it ============================ */ |
|
|
2083 |
|
|
2084 /* make sure degree list isn't empty */ |
|
|
2085 ASSERT (min_score >= 0) ; |
|
|
2086 ASSERT (min_score <= n_col) ; |
|
|
2087 ASSERT (head [min_score] >= EMPTY) ; |
|
|
2088 |
|
|
2089 #ifndef NDEBUG |
|
|
2090 for (debug_d = 0 ; debug_d < min_score ; debug_d++) |
|
|
2091 { |
|
|
2092 ASSERT (head [debug_d] == EMPTY) ; |
|
|
2093 } |
|
|
2094 #endif /* NDEBUG */ |
|
|
2095 |
|
|
2096 /* get pivot column from head of minimum degree list */ |
|
|
2097 while (head [min_score] == EMPTY && min_score < n_col) |
|
|
2098 { |
|
|
2099 min_score++ ; |
|
|
2100 } |
|
|
2101 pivot_col = head [min_score] ; |
|
|
2102 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; |
|
|
2103 next_col = Col [pivot_col].shared4.degree_next ; |
|
|
2104 head [min_score] = next_col ; |
|
|
2105 if (next_col != EMPTY) |
|
|
2106 { |
|
|
2107 Col [next_col].shared3.prev = EMPTY ; |
|
|
2108 } |
|
|
2109 |
|
|
2110 ASSERT (COL_IS_ALIVE (pivot_col)) ; |
|
|
2111 DEBUG3 (("Pivot col: %d\n", pivot_col)) ; |
|
|
2112 |
|
|
2113 /* remember score for defrag check */ |
|
|
2114 pivot_col_score = Col [pivot_col].shared2.score ; |
|
|
2115 |
|
|
2116 /* the pivot column is the kth column in the pivot order */ |
|
|
2117 Col [pivot_col].shared2.order = k ; |
|
|
2118 |
|
|
2119 /* increment order count by column thickness */ |
|
|
2120 pivot_col_thickness = Col [pivot_col].shared1.thickness ; |
|
|
2121 k += pivot_col_thickness ; |
|
|
2122 ASSERT (pivot_col_thickness > 0) ; |
|
|
2123 |
|
|
2124 /* === Garbage_collection, if necessary ============================= */ |
|
|
2125 |
|
|
2126 needed_memory = MIN (pivot_col_score, n_col - k) ; |
|
|
2127 if (pfree + needed_memory >= Alen) |
|
|
2128 { |
|
|
2129 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; |
|
|
2130 ngarbage++ ; |
|
|
2131 /* after garbage collection we will have enough */ |
|
|
2132 ASSERT (pfree + needed_memory < Alen) ; |
|
|
2133 /* garbage collection has wiped out the Row[].shared2.mark array */ |
|
|
2134 tag_mark = clear_mark (n_row, Row) ; |
|
|
2135 |
|
|
2136 #ifndef NDEBUG |
|
|
2137 debug_matrix (n_row, n_col, Row, Col, A) ; |
|
|
2138 #endif /* NDEBUG */ |
|
|
2139 } |
|
|
2140 |
|
|
2141 /* === Compute pivot row pattern ==================================== */ |
|
|
2142 |
|
|
2143 /* get starting location for this new merged row */ |
|
|
2144 pivot_row_start = pfree ; |
|
|
2145 |
|
|
2146 /* initialize new row counts to zero */ |
|
|
2147 pivot_row_degree = 0 ; |
|
|
2148 |
|
|
2149 /* tag pivot column as having been visited so it isn't included */ |
|
|
2150 /* in merged pivot row */ |
|
|
2151 Col [pivot_col].shared1.thickness = -pivot_col_thickness ; |
|
|
2152 |
|
|
2153 /* pivot row is the union of all rows in the pivot column pattern */ |
|
|
2154 cp = &A [Col [pivot_col].start] ; |
|
|
2155 cp_end = cp + Col [pivot_col].length ; |
|
|
2156 while (cp < cp_end) |
|
|
2157 { |
|
|
2158 /* get a row */ |
|
|
2159 row = *cp++ ; |
|
|
2160 DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; |
|
|
2161 /* skip if row is dead */ |
|
|
2162 if (ROW_IS_DEAD (row)) |
|
|
2163 { |
|
|
2164 continue ; |
|
|
2165 } |
|
|
2166 rp = &A [Row [row].start] ; |
|
|
2167 rp_end = rp + Row [row].length ; |
|
|
2168 while (rp < rp_end) |
|
|
2169 { |
|
|
2170 /* get a column */ |
|
|
2171 col = *rp++ ; |
|
|
2172 /* add the column, if alive and untagged */ |
|
|
2173 col_thickness = Col [col].shared1.thickness ; |
|
|
2174 if (col_thickness > 0 && COL_IS_ALIVE (col)) |
|
|
2175 { |
|
|
2176 /* tag column in pivot row */ |
|
|
2177 Col [col].shared1.thickness = -col_thickness ; |
|
|
2178 ASSERT (pfree < Alen) ; |
|
|
2179 /* place column in pivot row */ |
|
|
2180 A [pfree++] = col ; |
|
|
2181 pivot_row_degree += col_thickness ; |
|
|
2182 } |
|
|
2183 } |
|
|
2184 } |
|
|
2185 |
|
|
2186 /* clear tag on pivot column */ |
|
|
2187 Col [pivot_col].shared1.thickness = pivot_col_thickness ; |
|
|
2188 max_deg = MAX (max_deg, pivot_row_degree) ; |
|
|
2189 |
|
|
2190 #ifndef NDEBUG |
|
|
2191 DEBUG3 (("check2\n")) ; |
|
|
2192 debug_mark (n_row, Row, tag_mark, max_mark) ; |
|
|
2193 #endif /* NDEBUG */ |
|
|
2194 |
|
|
2195 /* === Kill all rows used to construct pivot row ==================== */ |
|
|
2196 |
|
|
2197 /* also kill pivot row, temporarily */ |
|
|
2198 cp = &A [Col [pivot_col].start] ; |
|
|
2199 cp_end = cp + Col [pivot_col].length ; |
|
|
2200 while (cp < cp_end) |
|
|
2201 { |
|
|
2202 /* may be killing an already dead row */ |
|
|
2203 row = *cp++ ; |
|
|
2204 DEBUG3 (("Kill row in pivot col: %d\n", row)) ; |
|
|
2205 KILL_ROW (row) ; |
|
|
2206 } |
|
|
2207 |
|
|
2208 /* === Select a row index to use as the new pivot row =============== */ |
|
|
2209 |
|
|
2210 pivot_row_length = pfree - pivot_row_start ; |
|
|
2211 if (pivot_row_length > 0) |
|
|
2212 { |
|
|
2213 /* pick the "pivot" row arbitrarily (first row in col) */ |
|
|
2214 pivot_row = A [Col [pivot_col].start] ; |
|
|
2215 DEBUG3 (("Pivotal row is %d\n", pivot_row)) ; |
|
|
2216 } |
|
|
2217 else |
|
|
2218 { |
|
|
2219 /* there is no pivot row, since it is of zero length */ |
|
|
2220 pivot_row = EMPTY ; |
|
|
2221 ASSERT (pivot_row_length == 0) ; |
|
|
2222 } |
|
|
2223 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; |
|
|
2224 |
|
|
2225 /* === Approximate degree computation =============================== */ |
|
|
2226 |
|
|
2227 /* Here begins the computation of the approximate degree. The column */ |
|
|
2228 /* score is the sum of the pivot row "length", plus the size of the */ |
|
|
2229 /* set differences of each row in the column minus the pattern of the */ |
|
|
2230 /* pivot row itself. The column ("thickness") itself is also */ |
|
|
2231 /* excluded from the column score (we thus use an approximate */ |
|
|
2232 /* external degree). */ |
|
|
2233 |
|
|
2234 /* The time taken by the following code (compute set differences, and */ |
|
|
2235 /* add them up) is proportional to the size of the data structure */ |
|
|
2236 /* being scanned - that is, the sum of the sizes of each column in */ |
|
|
2237 /* the pivot row. Thus, the amortized time to compute a column score */ |
|
|
2238 /* is proportional to the size of that column (where size, in this */ |
|
|
2239 /* context, is the column "length", or the number of row indices */ |
|
|
2240 /* in that column). The number of row indices in a column is */ |
|
|
2241 /* monotonically non-decreasing, from the length of the original */ |
|
|
2242 /* column on input to colamd. */ |
|
|
2243 |
|
|
2244 /* === Compute set differences ====================================== */ |
|
|
2245 |
|
|
2246 DEBUG3 (("** Computing set differences phase. **\n")) ; |
|
|
2247 |
|
|
2248 /* pivot row is currently dead - it will be revived later. */ |
|
|
2249 |
|
|
2250 DEBUG3 (("Pivot row: ")) ; |
|
|
2251 /* for each column in pivot row */ |
|
|
2252 rp = &A [pivot_row_start] ; |
|
|
2253 rp_end = rp + pivot_row_length ; |
|
|
2254 while (rp < rp_end) |
|
|
2255 { |
|
|
2256 col = *rp++ ; |
|
|
2257 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; |
|
|
2258 DEBUG3 (("Col: %d\n", col)) ; |
|
|
2259 |
|
|
2260 /* clear tags used to construct pivot row pattern */ |
|
|
2261 col_thickness = -Col [col].shared1.thickness ; |
|
|
2262 ASSERT (col_thickness > 0) ; |
|
|
2263 Col [col].shared1.thickness = col_thickness ; |
|
|
2264 |
|
|
2265 /* === Remove column from degree list =========================== */ |
|
|
2266 |
|
|
2267 cur_score = Col [col].shared2.score ; |
|
|
2268 prev_col = Col [col].shared3.prev ; |
|
|
2269 next_col = Col [col].shared4.degree_next ; |
|
|
2270 ASSERT (cur_score >= 0) ; |
|
|
2271 ASSERT (cur_score <= n_col) ; |
|
|
2272 ASSERT (cur_score >= EMPTY) ; |
|
|
2273 if (prev_col == EMPTY) |
|
|
2274 { |
|
|
2275 head [cur_score] = next_col ; |
|
|
2276 } |
|
|
2277 else |
|
|
2278 { |
|
|
2279 Col [prev_col].shared4.degree_next = next_col ; |
|
|
2280 } |
|
|
2281 if (next_col != EMPTY) |
|
|
2282 { |
|
|
2283 Col [next_col].shared3.prev = prev_col ; |
|
|
2284 } |
|
|
2285 |
|
|
2286 /* === Scan the column ========================================== */ |
|
|
2287 |
|
|
2288 cp = &A [Col [col].start] ; |
|
|
2289 cp_end = cp + Col [col].length ; |
|
|
2290 while (cp < cp_end) |
|
|
2291 { |
|
|
2292 /* get a row */ |
|
|
2293 row = *cp++ ; |
|
|
2294 row_mark = Row [row].shared2.mark ; |
|
|
2295 /* skip if dead */ |
|
|
2296 if (ROW_IS_MARKED_DEAD (row_mark)) |
|
|
2297 { |
|
|
2298 continue ; |
|
|
2299 } |
|
|
2300 ASSERT (row != pivot_row) ; |
|
|
2301 set_difference = row_mark - tag_mark ; |
|
|
2302 /* check if the row has been seen yet */ |
|
|
2303 if (set_difference < 0) |
|
|
2304 { |
|
|
2305 ASSERT (Row [row].shared1.degree <= max_deg) ; |
|
|
2306 set_difference = Row [row].shared1.degree ; |
|
|
2307 } |
|
|
2308 /* subtract column thickness from this row's set difference */ |
|
|
2309 set_difference -= col_thickness ; |
|
|
2310 ASSERT (set_difference >= 0) ; |
|
|
2311 /* absorb this row if the set difference becomes zero */ |
|
|
2312 if (set_difference == 0) |
|
|
2313 { |
|
|
2314 DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; |
|
|
2315 KILL_ROW (row) ; |
|
|
2316 } |
|
|
2317 else |
|
|
2318 { |
|
|
2319 /* save the new mark */ |
|
|
2320 Row [row].shared2.mark = set_difference + tag_mark ; |
|
|
2321 } |
|
|
2322 } |
|
|
2323 } |
|
|
2324 |
|
|
2325 #ifndef NDEBUG |
|
|
2326 debug_deg_lists (n_row, n_col, Row, Col, head, |
|
|
2327 min_score, n_col2-k-pivot_row_degree, max_deg) ; |
|
|
2328 #endif /* NDEBUG */ |
|
|
2329 |
|
|
2330 /* === Add up set differences for each column ======================= */ |
|
|
2331 |
|
|
2332 DEBUG3 (("** Adding set differences phase. **\n")) ; |
|
|
2333 |
|
|
2334 /* for each column in pivot row */ |
|
|
2335 rp = &A [pivot_row_start] ; |
|
|
2336 rp_end = rp + pivot_row_length ; |
|
|
2337 while (rp < rp_end) |
|
|
2338 { |
|
|
2339 /* get a column */ |
|
|
2340 col = *rp++ ; |
|
|
2341 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; |
|
|
2342 hash = 0 ; |
|
|
2343 cur_score = 0 ; |
|
|
2344 cp = &A [Col [col].start] ; |
|
|
2345 /* compact the column */ |
|
|
2346 new_cp = cp ; |
|
|
2347 cp_end = cp + Col [col].length ; |
|
|
2348 |
|
|
2349 DEBUG4 (("Adding set diffs for Col: %d.\n", col)) ; |
|
|
2350 |
|
|
2351 while (cp < cp_end) |
|
|
2352 { |
|
|
2353 /* get a row */ |
|
|
2354 row = *cp++ ; |
|
|
2355 ASSERT(row >= 0 && row < n_row) ; |
|
|
2356 row_mark = Row [row].shared2.mark ; |
|
|
2357 /* skip if dead */ |
|
|
2358 if (ROW_IS_MARKED_DEAD (row_mark)) |
|
|
2359 { |
|
|
2360 continue ; |
|
|
2361 } |
|
|
2362 ASSERT (row_mark > tag_mark) ; |
|
|
2363 /* compact the column */ |
|
|
2364 *new_cp++ = row ; |
|
|
2365 /* compute hash function */ |
|
|
2366 hash += row ; |
|
|
2367 /* add set difference */ |
|
|
2368 cur_score += row_mark - tag_mark ; |
|
|
2369 /* integer overflow... */ |
|
|
2370 cur_score = MIN (cur_score, n_col) ; |
|
|
2371 } |
|
|
2372 |
|
|
2373 /* recompute the column's length */ |
|
|
2374 Col [col].length = (int) (new_cp - &A [Col [col].start]) ; |
|
|
2375 |
|
|
2376 /* === Further mass elimination ================================= */ |
|
|
2377 |
|
|
2378 if (Col [col].length == 0) |
|
|
2379 { |
|
|
2380 DEBUG4 (("further mass elimination. Col: %d\n", col)) ; |
|
|
2381 /* nothing left but the pivot row in this column */ |
|
|
2382 KILL_PRINCIPAL_COL (col) ; |
|
|
2383 pivot_row_degree -= Col [col].shared1.thickness ; |
|
|
2384 ASSERT (pivot_row_degree >= 0) ; |
|
|
2385 /* order it */ |
|
|
2386 Col [col].shared2.order = k ; |
|
|
2387 /* increment order count by column thickness */ |
|
|
2388 k += Col [col].shared1.thickness ; |
|
|
2389 } |
|
|
2390 else |
|
|
2391 { |
|
|
2392 /* === Prepare for supercolumn detection ==================== */ |
|
|
2393 |
|
|
2394 DEBUG4 (("Preparing supercol detection for Col: %d.\n", col)) ; |
|
|
2395 |
|
|
2396 /* save score so far */ |
|
|
2397 Col [col].shared2.score = cur_score ; |
|
|
2398 |
|
|
2399 /* add column to hash table, for supercolumn detection */ |
|
|
2400 hash %= n_col + 1 ; |
|
|
2401 |
|
|
2402 DEBUG4 ((" Hash = %d, n_col = %d.\n", hash, n_col)) ; |
|
|
2403 ASSERT (hash <= n_col) ; |
|
|
2404 |
|
|
2405 head_column = head [hash] ; |
|
|
2406 if (head_column > EMPTY) |
|
|
2407 { |
|
|
2408 /* degree list "hash" is non-empty, use prev (shared3) of */ |
|
|
2409 /* first column in degree list as head of hash bucket */ |
|
|
2410 first_col = Col [head_column].shared3.headhash ; |
|
|
2411 Col [head_column].shared3.headhash = col ; |
|
|
2412 } |
|
|
2413 else |
|
|
2414 { |
|
|
2415 /* degree list "hash" is empty, use head as hash bucket */ |
|
|
2416 first_col = - (head_column + 2) ; |
|
|
2417 head [hash] = - (col + 2) ; |
|
|
2418 } |
|
|
2419 Col [col].shared4.hash_next = first_col ; |
|
|
2420 |
|
|
2421 /* save hash function in Col [col].shared3.hash */ |
|
|
2422 Col [col].shared3.hash = (int) hash ; |
|
|
2423 ASSERT (COL_IS_ALIVE (col)) ; |
|
|
2424 } |
|
|
2425 } |
|
|
2426 |
|
|
2427 /* The approximate external column degree is now computed. */ |
|
|
2428 |
|
|
2429 /* === Supercolumn detection ======================================== */ |
|
|
2430 |
|
|
2431 DEBUG3 (("** Supercolumn detection phase. **\n")) ; |
|
|
2432 |
|
|
2433 detect_super_cols ( |
|
|
2434 |
|
|
2435 #ifndef NDEBUG |
|
|
2436 n_col, Row, |
|
|
2437 #endif /* NDEBUG */ |
|
|
2438 |
|
|
2439 Col, A, head, pivot_row_start, pivot_row_length) ; |
|
|
2440 |
|
|
2441 /* === Kill the pivotal column ====================================== */ |
|
|
2442 |
|
|
2443 KILL_PRINCIPAL_COL (pivot_col) ; |
|
|
2444 |
|
|
2445 /* === Clear mark =================================================== */ |
|
|
2446 |
|
|
2447 tag_mark += (max_deg + 1) ; |
|
|
2448 if (tag_mark >= max_mark) |
|
|
2449 { |
|
|
2450 DEBUG2 (("clearing tag_mark\n")) ; |
|
|
2451 tag_mark = clear_mark (n_row, Row) ; |
|
|
2452 } |
|
|
2453 |
|
|
2454 #ifndef NDEBUG |
|
|
2455 DEBUG3 (("check3\n")) ; |
|
|
2456 debug_mark (n_row, Row, tag_mark, max_mark) ; |
|
|
2457 #endif /* NDEBUG */ |
|
|
2458 |
|
|
2459 /* === Finalize the new pivot row, and column scores ================ */ |
|
|
2460 |
|
|
2461 DEBUG3 (("** Finalize scores phase. **\n")) ; |
|
|
2462 |
|
|
2463 /* for each column in pivot row */ |
|
|
2464 rp = &A [pivot_row_start] ; |
|
|
2465 /* compact the pivot row */ |
|
|
2466 new_rp = rp ; |
|
|
2467 rp_end = rp + pivot_row_length ; |
|
|
2468 while (rp < rp_end) |
|
|
2469 { |
|
|
2470 col = *rp++ ; |
|
|
2471 /* skip dead columns */ |
|
|
2472 if (COL_IS_DEAD (col)) |
|
|
2473 { |
|
|
2474 continue ; |
|
|
2475 } |
|
|
2476 *new_rp++ = col ; |
|
|
2477 /* add new pivot row to column */ |
|
|
2478 A [Col [col].start + (Col [col].length++)] = pivot_row ; |
|
|
2479 |
|
|
2480 /* retrieve score so far and add on pivot row's degree. */ |
|
|
2481 /* (we wait until here for this in case the pivot */ |
|
|
2482 /* row's degree was reduced due to mass elimination). */ |
|
|
2483 cur_score = Col [col].shared2.score + pivot_row_degree ; |
|
|
2484 |
|
|
2485 /* calculate the max possible score as the number of */ |
|
|
2486 /* external columns minus the 'k' value minus the */ |
|
|
2487 /* columns thickness */ |
|
|
2488 max_score = n_col - k - Col [col].shared1.thickness ; |
|
|
2489 |
|
|
2490 /* make the score the external degree of the union-of-rows */ |
|
|
2491 cur_score -= Col [col].shared1.thickness ; |
|
|
2492 |
|
|
2493 /* make sure score is less or equal than the max score */ |
|
|
2494 cur_score = MIN (cur_score, max_score) ; |
|
|
2495 ASSERT (cur_score >= 0) ; |
|
|
2496 |
|
|
2497 /* store updated score */ |
|
|
2498 Col [col].shared2.score = cur_score ; |
|
|
2499 |
|
|
2500 /* === Place column back in degree list ========================= */ |
|
|
2501 |
|
|
2502 ASSERT (min_score >= 0) ; |
|
|
2503 ASSERT (min_score <= n_col) ; |
|
|
2504 ASSERT (cur_score >= 0) ; |
|
|
2505 ASSERT (cur_score <= n_col) ; |
|
|
2506 ASSERT (head [cur_score] >= EMPTY) ; |
|
|
2507 next_col = head [cur_score] ; |
|
|
2508 Col [col].shared4.degree_next = next_col ; |
|
|
2509 Col [col].shared3.prev = EMPTY ; |
|
|
2510 if (next_col != EMPTY) |
|
|
2511 { |
|
|
2512 Col [next_col].shared3.prev = col ; |
|
|
2513 } |
|
|
2514 head [cur_score] = col ; |
|
|
2515 |
|
|
2516 /* see if this score is less than current min */ |
|
|
2517 min_score = MIN (min_score, cur_score) ; |
|
|
2518 |
|
|
2519 } |
|
|
2520 |
|
|
2521 #ifndef NDEBUG |
|
|
2522 debug_deg_lists (n_row, n_col, Row, Col, head, |
|
|
2523 min_score, n_col2-k, max_deg) ; |
|
|
2524 #endif /* NDEBUG */ |
|
|
2525 |
|
|
2526 /* === Resurrect the new pivot row ================================== */ |
|
|
2527 |
|
|
2528 if (pivot_row_degree > 0) |
|
|
2529 { |
|
|
2530 /* update pivot row length to reflect any cols that were killed */ |
|
|
2531 /* during super-col detection and mass elimination */ |
|
|
2532 Row [pivot_row].start = pivot_row_start ; |
|
|
2533 Row [pivot_row].length = (int) (new_rp - &A[pivot_row_start]) ; |
|
|
2534 Row [pivot_row].shared1.degree = pivot_row_degree ; |
|
|
2535 Row [pivot_row].shared2.mark = 0 ; |
|
|
2536 /* pivot row is no longer dead */ |
|
|
2537 } |
|
|
2538 } |
|
|
2539 |
|
|
2540 /* === All principal columns have now been ordered ====================== */ |
|
|
2541 |
|
|
2542 return (ngarbage) ; |
|
|
2543 } |
|
|
2544 |
|
|
2545 |
|
|
2546 /* ========================================================================== */ |
|
|
2547 /* === order_children ======================================================= */ |
|
|
2548 /* ========================================================================== */ |
|
|
2549 |
|
|
2550 /* |
|
|
2551 The find_ordering routine has ordered all of the principal columns (the |
|
|
2552 representatives of the supercolumns). The non-principal columns have not |
|
|
2553 yet been ordered. This routine orders those columns by walking up the |
|
|
2554 parent tree (a column is a child of the column which absorbed it). The |
|
|
2555 final permutation vector is then placed in p [0 ... n_col-1], with p [0] |
|
|
2556 being the first column, and p [n_col-1] being the last. It doesn't look |
|
|
2557 like it at first glance, but be assured that this routine takes time linear |
|
|
2558 in the number of columns. Although not immediately obvious, the time |
|
|
2559 taken by this routine is O (n_col), that is, linear in the number of |
|
|
2560 columns. Not user-callable. |
|
|
2561 */ |
|
|
2562 |
|
|
2563 PRIVATE void order_children |
|
|
2564 ( |
|
|
2565 /* === Parameters ======================================================= */ |
|
|
2566 |
|
|
2567 int n_col, /* number of columns of A */ |
|
|
2568 Colamd_Col Col [], /* of size n_col+1 */ |
|
|
2569 int p [] /* p [0 ... n_col-1] is the column permutation*/ |
|
|
2570 ) |
|
|
2571 { |
|
|
2572 /* === Local variables ================================================== */ |
|
|
2573 |
|
|
2574 int i ; /* loop counter for all columns */ |
|
|
2575 int c ; /* column index */ |
|
|
2576 int parent ; /* index of column's parent */ |
|
|
2577 int order ; /* column's order */ |
|
|
2578 |
|
|
2579 /* === Order each non-principal column ================================== */ |
|
|
2580 |
|
|
2581 for (i = 0 ; i < n_col ; i++) |
|
|
2582 { |
|
|
2583 /* find an un-ordered non-principal column */ |
|
|
2584 ASSERT (COL_IS_DEAD (i)) ; |
|
|
2585 if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EMPTY) |
|
|
2586 { |
|
|
2587 parent = i ; |
|
|
2588 /* once found, find its principal parent */ |
|
|
2589 do |
|
|
2590 { |
|
|
2591 parent = Col [parent].shared1.parent ; |
|
|
2592 } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; |
|
|
2593 |
|
|
2594 /* now, order all un-ordered non-principal columns along path */ |
|
|
2595 /* to this parent. collapse tree at the same time */ |
|
|
2596 c = i ; |
|
|
2597 /* get order of parent */ |
|
|
2598 order = Col [parent].shared2.order ; |
|
|
2599 |
|
|
2600 do |
|
|
2601 { |
|
|
2602 ASSERT (Col [c].shared2.order == EMPTY) ; |
|
|
2603 |
|
|
2604 /* order this column */ |
|
|
2605 Col [c].shared2.order = order++ ; |
|
|
2606 /* collaps tree */ |
|
|
2607 Col [c].shared1.parent = parent ; |
|
|
2608 |
|
|
2609 /* get immediate parent of this column */ |
|
|
2610 c = Col [c].shared1.parent ; |
|
|
2611 |
|
|
2612 /* continue until we hit an ordered column. There are */ |
|
|
2613 /* guarranteed not to be anymore unordered columns */ |
|
|
2614 /* above an ordered column */ |
|
|
2615 } while (Col [c].shared2.order == EMPTY) ; |
|
|
2616 |
|
|
2617 /* re-order the super_col parent to largest order for this group */ |
|
|
2618 Col [parent].shared2.order = order ; |
|
|
2619 } |
|
|
2620 } |
|
|
2621 |
|
|
2622 /* === Generate the permutation ========================================= */ |
|
|
2623 |
|
|
2624 for (c = 0 ; c < n_col ; c++) |
|
|
2625 { |
|
|
2626 p [Col [c].shared2.order] = c ; |
|
|
2627 } |
|
|
2628 } |
|
|
2629 |
|
|
2630 |
|
|
2631 /* ========================================================================== */ |
|
|
2632 /* === detect_super_cols ==================================================== */ |
|
|
2633 /* ========================================================================== */ |
|
|
2634 |
|
|
2635 /* |
|
|
2636 Detects supercolumns by finding matches between columns in the hash buckets. |
|
|
2637 Check amongst columns in the set A [row_start ... row_start + row_length-1]. |
|
|
2638 The columns under consideration are currently *not* in the degree lists, |
|
|
2639 and have already been placed in the hash buckets. |
|
|
2640 |
|
|
2641 The hash bucket for columns whose hash function is equal to h is stored |
|
|
2642 as follows: |
|
|
2643 |
|
|
2644 if head [h] is >= 0, then head [h] contains a degree list, so: |
|
|
2645 |
|
|
2646 head [h] is the first column in degree bucket h. |
|
|
2647 Col [head [h]].headhash gives the first column in hash bucket h. |
|
|
2648 |
|
|
2649 otherwise, the degree list is empty, and: |
|
|
2650 |
|
|
2651 -(head [h] + 2) is the first column in hash bucket h. |
|
|
2652 |
|
|
2653 For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous |
|
|
2654 column" pointer. Col [c].shared3.hash is used instead as the hash number |
|
|
2655 for that column. The value of Col [c].shared4.hash_next is the next column |
|
|
2656 in the same hash bucket. |
|
|
2657 |
|
|
2658 Assuming no, or "few" hash collisions, the time taken by this routine is |
|
|
2659 linear in the sum of the sizes (lengths) of each column whose score has |
|
|
2660 just been computed in the approximate degree computation. |
|
|
2661 Not user-callable. |
|
|
2662 */ |
|
|
2663 |
|
|
2664 PRIVATE void detect_super_cols |
|
|
2665 ( |
|
|
2666 /* === Parameters ======================================================= */ |
|
|
2667 |
|
|
2668 #ifndef NDEBUG |
|
|
2669 /* these two parameters are only needed when debugging is enabled: */ |
|
|
2670 int n_col, /* number of columns of A */ |
|
|
2671 Colamd_Row Row [], /* of size n_row+1 */ |
|
|
2672 #endif /* NDEBUG */ |
|
|
2673 |
|
|
2674 Colamd_Col Col [], /* of size n_col+1 */ |
|
|
2675 int A [], /* row indices of A */ |
|
|
2676 int head [], /* head of degree lists and hash buckets */ |
|
|
2677 int row_start, /* pointer to set of columns to check */ |
|
|
2678 int row_length /* number of columns to check */ |
|
|
2679 ) |
|
|
2680 { |
|
|
2681 /* === Local variables ================================================== */ |
|
|
2682 |
|
|
2683 int hash ; /* hash value for a column */ |
|
|
2684 int *rp ; /* pointer to a row */ |
|
|
2685 int c ; /* a column index */ |
|
|
2686 int super_c ; /* column index of the column to absorb into */ |
|
|
2687 int *cp1 ; /* column pointer for column super_c */ |
|
|
2688 int *cp2 ; /* column pointer for column c */ |
|
|
2689 int length ; /* length of column super_c */ |
|
|
2690 int prev_c ; /* column preceding c in hash bucket */ |
|
|
2691 int i ; /* loop counter */ |
|
|
2692 int *rp_end ; /* pointer to the end of the row */ |
|
|
2693 int col ; /* a column index in the row to check */ |
|
|
2694 int head_column ; /* first column in hash bucket or degree list */ |
|
|
2695 int first_col ; /* first column in hash bucket */ |
|
|
2696 |
|
|
2697 /* === Consider each column in the row ================================== */ |
|
|
2698 |
|
|
2699 rp = &A [row_start] ; |
|
|
2700 rp_end = rp + row_length ; |
|
|
2701 while (rp < rp_end) |
|
|
2702 { |
|
|
2703 col = *rp++ ; |
|
|
2704 if (COL_IS_DEAD (col)) |
|
|
2705 { |
|
|
2706 continue ; |
|
|
2707 } |
|
|
2708 |
|
|
2709 /* get hash number for this column */ |
|
|
2710 hash = Col [col].shared3.hash ; |
|
|
2711 ASSERT (hash <= n_col) ; |
|
|
2712 |
|
|
2713 /* === Get the first column in this hash bucket ===================== */ |
|
|
2714 |
|
|
2715 head_column = head [hash] ; |
|
|
2716 if (head_column > EMPTY) |
|
|
2717 { |
|
|
2718 first_col = Col [head_column].shared3.headhash ; |
|
|
2719 } |
|
|
2720 else |
|
|
2721 { |
|
|
2722 first_col = - (head_column + 2) ; |
|
|
2723 } |
|
|
2724 |
|
|
2725 /* === Consider each column in the hash bucket ====================== */ |
|
|
2726 |
|
|
2727 for (super_c = first_col ; super_c != EMPTY ; |
|
|
2728 super_c = Col [super_c].shared4.hash_next) |
|
|
2729 { |
|
|
2730 ASSERT (COL_IS_ALIVE (super_c)) ; |
|
|
2731 ASSERT (Col [super_c].shared3.hash == hash) ; |
|
|
2732 length = Col [super_c].length ; |
|
|
2733 |
|
|
2734 /* prev_c is the column preceding column c in the hash bucket */ |
|
|
2735 prev_c = super_c ; |
|
|
2736 |
|
|
2737 /* === Compare super_c with all columns after it ================ */ |
|
|
2738 |
|
|
2739 for (c = Col [super_c].shared4.hash_next ; |
|
|
2740 c != EMPTY ; c = Col [c].shared4.hash_next) |
|
|
2741 { |
|
|
2742 ASSERT (c != super_c) ; |
|
|
2743 ASSERT (COL_IS_ALIVE (c)) ; |
|
|
2744 ASSERT (Col [c].shared3.hash == hash) ; |
|
|
2745 |
|
|
2746 /* not identical if lengths or scores are different */ |
|
|
2747 if (Col [c].length != length || |
|
|
2748 Col [c].shared2.score != Col [super_c].shared2.score) |
|
|
2749 { |
|
|
2750 prev_c = c ; |
|
|
2751 continue ; |
|
|
2752 } |
|
|
2753 |
|
|
2754 /* compare the two columns */ |
|
|
2755 cp1 = &A [Col [super_c].start] ; |
|
|
2756 cp2 = &A [Col [c].start] ; |
|
|
2757 |
|
|
2758 for (i = 0 ; i < length ; i++) |
|
|
2759 { |
|
|
2760 /* the columns are "clean" (no dead rows) */ |
|
|
2761 ASSERT (ROW_IS_ALIVE (*cp1)) ; |
|
|
2762 ASSERT (ROW_IS_ALIVE (*cp2)) ; |
|
|
2763 /* row indices will same order for both supercols, */ |
|
|
2764 /* no gather scatter nessasary */ |
|
|
2765 if (*cp1++ != *cp2++) |
|
|
2766 { |
|
|
2767 break ; |
|
|
2768 } |
|
|
2769 } |
|
|
2770 |
|
|
2771 /* the two columns are different if the for-loop "broke" */ |
|
|
2772 if (i != length) |
|
|
2773 { |
|
|
2774 prev_c = c ; |
|
|
2775 continue ; |
|
|
2776 } |
|
|
2777 |
|
|
2778 /* === Got it! two columns are identical =================== */ |
|
|
2779 |
|
|
2780 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ; |
|
|
2781 |
|
|
2782 Col [super_c].shared1.thickness += Col [c].shared1.thickness ; |
|
|
2783 Col [c].shared1.parent = super_c ; |
|
|
2784 KILL_NON_PRINCIPAL_COL (c) ; |
|
|
2785 /* order c later, in order_children() */ |
|
|
2786 Col [c].shared2.order = EMPTY ; |
|
|
2787 /* remove c from hash bucket */ |
|
|
2788 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; |
|
|
2789 } |
|
|
2790 } |
|
|
2791 |
|
|
2792 /* === Empty this hash bucket ======================================= */ |
|
|
2793 |
|
|
2794 if (head_column > EMPTY) |
|
|
2795 { |
|
|
2796 /* corresponding degree list "hash" is not empty */ |
|
|
2797 Col [head_column].shared3.headhash = EMPTY ; |
|
|
2798 } |
|
|
2799 else |
|
|
2800 { |
|
|
2801 /* corresponding degree list "hash" is empty */ |
|
|
2802 head [hash] = EMPTY ; |
|
|
2803 } |
|
|
2804 } |
|
|
2805 } |
|
|
2806 |
|
|
2807 |
|
|
2808 /* ========================================================================== */ |
|
|
2809 /* === garbage_collection =================================================== */ |
|
|
2810 /* ========================================================================== */ |
|
|
2811 |
|
|
2812 /* |
|
|
2813 Defragments and compacts columns and rows in the workspace A. Used when |
|
|
2814 all avaliable memory has been used while performing row merging. Returns |
|
|
2815 the index of the first free position in A, after garbage collection. The |
|
|
2816 time taken by this routine is linear is the size of the array A, which is |
|
|
2817 itself linear in the number of nonzeros in the input matrix. |
|
|
2818 Not user-callable. |
|
|
2819 */ |
|
|
2820 |
|
|
2821 PRIVATE int garbage_collection /* returns the new value of pfree */ |
|
|
2822 ( |
|
|
2823 /* === Parameters ======================================================= */ |
|
|
2824 |
|
|
2825 int n_row, /* number of rows */ |
|
|
2826 int n_col, /* number of columns */ |
|
|
2827 Colamd_Row Row [], /* row info */ |
|
|
2828 Colamd_Col Col [], /* column info */ |
|
|
2829 int A [], /* A [0 ... Alen-1] holds the matrix */ |
|
|
2830 int *pfree /* &A [0] ... pfree is in use */ |
|
|
2831 ) |
|
|
2832 { |
|
|
2833 /* === Local variables ================================================== */ |
|
|
2834 |
|
|
2835 int *psrc ; /* source pointer */ |
|
|
2836 int *pdest ; /* destination pointer */ |
|
|
2837 int j ; /* counter */ |
|
|
2838 int r ; /* a row index */ |
|
|
2839 int c ; /* a column index */ |
|
|
2840 int length ; /* length of a row or column */ |
|
|
2841 |
|
|
2842 #ifndef NDEBUG |
|
|
2843 int debug_rows ; |
|
|
2844 DEBUG2 (("Defrag..\n")) ; |
|
|
2845 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ; |
|
|
2846 debug_rows = 0 ; |
|
|
2847 #endif /* NDEBUG */ |
|
|
2848 |
|
|
2849 /* === Defragment the columns =========================================== */ |
|
|
2850 |
|
|
2851 pdest = &A[0] ; |
|
|
2852 for (c = 0 ; c < n_col ; c++) |
|
|
2853 { |
|
|
2854 if (COL_IS_ALIVE (c)) |
|
|
2855 { |
|
|
2856 psrc = &A [Col [c].start] ; |
|
|
2857 |
|
|
2858 /* move and compact the column */ |
|
|
2859 ASSERT (pdest <= psrc) ; |
|
|
2860 Col [c].start = (int) (pdest - &A [0]) ; |
|
|
2861 length = Col [c].length ; |
|
|
2862 for (j = 0 ; j < length ; j++) |
|
|
2863 { |
|
|
2864 r = *psrc++ ; |
|
|
2865 if (ROW_IS_ALIVE (r)) |
|
|
2866 { |
|
|
2867 *pdest++ = r ; |
|
|
2868 } |
|
|
2869 } |
|
|
2870 Col [c].length = (int) (pdest - &A [Col [c].start]) ; |
|
|
2871 } |
|
|
2872 } |
|
|
2873 |
|
|
2874 /* === Prepare to defragment the rows =================================== */ |
|
|
2875 |
|
|
2876 for (r = 0 ; r < n_row ; r++) |
|
|
2877 { |
|
|
2878 if (ROW_IS_ALIVE (r)) |
|
|
2879 { |
|
|
2880 if (Row [r].length == 0) |
|
|
2881 { |
|
|
2882 /* this row is of zero length. cannot compact it, so kill it */ |
|
|
2883 DEBUG3 (("Defrag row kill\n")) ; |
|
|
2884 KILL_ROW (r) ; |
|
|
2885 } |
|
|
2886 else |
|
|
2887 { |
|
|
2888 /* save first column index in Row [r].shared2.first_column */ |
|
|
2889 psrc = &A [Row [r].start] ; |
|
|
2890 Row [r].shared2.first_column = *psrc ; |
|
|
2891 ASSERT (ROW_IS_ALIVE (r)) ; |
|
|
2892 /* flag the start of the row with the one's complement of row */ |
|
|
2893 *psrc = ONES_COMPLEMENT (r) ; |
|
|
2894 |
|
|
2895 #ifndef NDEBUG |
|
|
2896 debug_rows++ ; |
|
|
2897 #endif /* NDEBUG */ |
|
|
2898 |
|
|
2899 } |
|
|
2900 } |
|
|
2901 } |
|
|
2902 |
|
|
2903 /* === Defragment the rows ============================================== */ |
|
|
2904 |
|
|
2905 psrc = pdest ; |
|
|
2906 while (psrc < pfree) |
|
|
2907 { |
|
|
2908 /* find a negative number ... the start of a row */ |
|
|
2909 if (*psrc++ < 0) |
|
|
2910 { |
|
|
2911 psrc-- ; |
|
|
2912 /* get the row index */ |
|
|
2913 r = ONES_COMPLEMENT (*psrc) ; |
|
|
2914 ASSERT (r >= 0 && r < n_row) ; |
|
|
2915 /* restore first column index */ |
|
|
2916 *psrc = Row [r].shared2.first_column ; |
|
|
2917 ASSERT (ROW_IS_ALIVE (r)) ; |
|
|
2918 |
|
|
2919 /* move and compact the row */ |
|
|
2920 ASSERT (pdest <= psrc) ; |
|
|
2921 Row [r].start = (int) (pdest - &A [0]) ; |
|
|
2922 length = Row [r].length ; |
|
|
2923 for (j = 0 ; j < length ; j++) |
|
|
2924 { |
|
|
2925 c = *psrc++ ; |
|
|
2926 if (COL_IS_ALIVE (c)) |
|
|
2927 { |
|
|
2928 *pdest++ = c ; |
|
|
2929 } |
|
|
2930 } |
|
|
2931 Row [r].length = (int) (pdest - &A [Row [r].start]) ; |
|
|
2932 |
|
|
2933 #ifndef NDEBUG |
|
|
2934 debug_rows-- ; |
|
|
2935 #endif /* NDEBUG */ |
|
|
2936 |
|
|
2937 } |
|
|
2938 } |
|
|
2939 /* ensure we found all the rows */ |
|
|
2940 ASSERT (debug_rows == 0) ; |
|
|
2941 |
|
|
2942 /* === Return the new value of pfree ==================================== */ |
|
|
2943 |
|
|
2944 return ((int) (pdest - &A [0])) ; |
|
|
2945 } |
|
|
2946 |
|
|
2947 |
|
|
2948 /* ========================================================================== */ |
|
|
2949 /* === clear_mark =========================================================== */ |
|
|
2950 /* ========================================================================== */ |
|
|
2951 |
|
|
2952 /* |
|
|
2953 Clears the Row [].shared2.mark array, and returns the new tag_mark. |
|
|
2954 Return value is the new tag_mark. Not user-callable. |
|
|
2955 */ |
|
|
2956 |
|
|
2957 PRIVATE int clear_mark /* return the new value for tag_mark */ |
|
|
2958 ( |
|
|
2959 /* === Parameters ======================================================= */ |
|
|
2960 |
|
|
2961 int n_row, /* number of rows in A */ |
|
|
2962 Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ |
|
|
2963 ) |
|
|
2964 { |
|
|
2965 /* === Local variables ================================================== */ |
|
|
2966 |
|
|
2967 int r ; |
|
|
2968 |
|
|
2969 for (r = 0 ; r < n_row ; r++) |
|
|
2970 { |
|
|
2971 if (ROW_IS_ALIVE (r)) |
|
|
2972 { |
|
|
2973 Row [r].shared2.mark = 0 ; |
|
|
2974 } |
|
|
2975 } |
|
|
2976 return (1) ; |
|
|
2977 } |
|
|
2978 |
|
|
2979 |
|
|
2980 /* ========================================================================== */ |
|
|
2981 /* === print_report ========================================================= */ |
|
|
2982 /* ========================================================================== */ |
|
|
2983 |
|
|
2984 PRIVATE void print_report |
|
|
2985 ( |
|
|
2986 char *method, |
|
|
2987 int stats [COLAMD_STATS] |
|
|
2988 ) |
|
|
2989 { |
|
|
2990 |
|
|
2991 int i1, i2, i3 ; |
|
|
2992 |
|
|
2993 if (!stats) |
|
|
2994 { |
|
|
2995 PRINTF ("%s: No statistics available.\n", method) ; |
|
|
2996 return ; |
|
|
2997 } |
|
|
2998 |
|
|
2999 i1 = stats [COLAMD_INFO1] ; |
|
|
3000 i2 = stats [COLAMD_INFO2] ; |
|
|
3001 i3 = stats [COLAMD_INFO3] ; |
|
|
3002 |
|
|
3003 if (stats [COLAMD_STATUS] >= 0) |
|
|
3004 { |
|
|
3005 PRINTF ("%s: OK. ", method) ; |
|
|
3006 } |
|
|
3007 else |
|
|
3008 { |
|
|
3009 PRINTF ("%s: ERROR. ", method) ; |
|
|
3010 } |
|
|
3011 |
|
|
3012 switch (stats [COLAMD_STATUS]) |
|
|
3013 { |
|
|
3014 |
|
|
3015 case COLAMD_OK_BUT_JUMBLED: |
|
|
3016 |
|
|
3017 PRINTF ("Matrix has unsorted or duplicate row indices.\n") ; |
|
|
3018 |
|
|
3019 PRINTF ("%s: number of duplicate or out-of-order row indices: %d\n", |
|
|
3020 method, i3) ; |
|
|
3021 |
|
|
3022 PRINTF ("%s: last seen duplicate or out-of-order row index: %d\n", |
|
|
3023 method, INDEX (i2)) ; |
|
|
3024 |
|
|
3025 PRINTF ("%s: last seen in column: %d", |
|
|
3026 method, INDEX (i1)) ; |
|
|
3027 |
|
|
3028 /* no break - fall through to next case instead */ |
|
|
3029 |
|
|
3030 case COLAMD_OK: |
|
|
3031 |
|
|
3032 PRINTF ("\n") ; |
|
|
3033 |
|
|
3034 PRINTF ("%s: number of dense or empty rows ignored: %d\n", |
|
|
3035 method, stats [COLAMD_DENSE_ROW]) ; |
|
|
3036 |
|
|
3037 PRINTF ("%s: number of dense or empty columns ignored: %d\n", |
|
|
3038 method, stats [COLAMD_DENSE_COL]) ; |
|
|
3039 |
|
|
3040 PRINTF ("%s: number of garbage collections performed: %d\n", |
|
|
3041 method, stats [COLAMD_DEFRAG_COUNT]) ; |
|
|
3042 break ; |
|
|
3043 |
|
|
3044 case COLAMD_ERROR_A_not_present: |
|
|
3045 |
|
|
3046 PRINTF ("Array A (row indices of matrix) not present.\n") ; |
|
|
3047 break ; |
|
|
3048 |
|
|
3049 case COLAMD_ERROR_p_not_present: |
|
|
3050 |
|
|
3051 PRINTF ("Array p (column pointers for matrix) not present.\n") ; |
|
|
3052 break ; |
|
|
3053 |
|
|
3054 case COLAMD_ERROR_nrow_negative: |
|
|
3055 |
|
|
3056 PRINTF ("Invalid number of rows (%d).\n", i1) ; |
|
|
3057 break ; |
|
|
3058 |
|
|
3059 case COLAMD_ERROR_ncol_negative: |
|
|
3060 |
|
|
3061 PRINTF ("Invalid number of columns (%d).\n", i1) ; |
|
|
3062 break ; |
|
|
3063 |
|
|
3064 case COLAMD_ERROR_nnz_negative: |
|
|
3065 |
|
|
3066 PRINTF ("Invalid number of nonzero entries (%d).\n", i1) ; |
|
|
3067 break ; |
|
|
3068 |
|
|
3069 case COLAMD_ERROR_p0_nonzero: |
|
|
3070 |
|
|
3071 PRINTF ("Invalid column pointer, p [0] = %d, must be zero.\n", i1) ; |
|
|
3072 break ; |
|
|
3073 |
|
|
3074 case COLAMD_ERROR_A_too_small: |
|
|
3075 |
|
|
3076 PRINTF ("Array A too small.\n") ; |
|
|
3077 PRINTF (" Need Alen >= %d, but given only Alen = %d.\n", |
|
|
3078 i1, i2) ; |
|
|
3079 break ; |
|
|
3080 |
|
|
3081 case COLAMD_ERROR_col_length_negative: |
|
|
3082 |
|
|
3083 PRINTF |
|
|
3084 ("Column %d has a negative number of nonzero entries (%d).\n", |
|
|
3085 INDEX (i1), i2) ; |
|
|
3086 break ; |
|
|
3087 |
|
|
3088 case COLAMD_ERROR_row_index_out_of_bounds: |
|
|
3089 |
|
|
3090 PRINTF |
|
|
3091 ("Row index (row %d) out of bounds (%d to %d) in column %d.\n", |
|
|
3092 INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1)) ; |
|
|
3093 break ; |
|
|
3094 |
|
|
3095 case COLAMD_ERROR_out_of_memory: |
|
|
3096 |
|
|
3097 PRINTF ("Out of memory.\n") ; |
|
|
3098 break ; |
|
|
3099 |
|
|
3100 case COLAMD_ERROR_internal_error: |
|
|
3101 |
|
|
3102 /* if this happens, there is a bug in the code */ |
|
|
3103 PRINTF |
|
|
3104 ("Internal error! Please contact authors (davis@cise.ufl.edu).\n") ; |
|
|
3105 break ; |
|
|
3106 } |
|
|
3107 } |
|
|
3108 |
|
|
3109 |
|
|
3110 |
|
|
3111 |
|
|
3112 /* ========================================================================== */ |
|
|
3113 /* === colamd debugging routines ============================================ */ |
|
|
3114 /* ========================================================================== */ |
|
|
3115 |
|
|
3116 /* When debugging is disabled, the remainder of this file is ignored. */ |
|
|
3117 |
|
|
3118 #ifndef NDEBUG |
|
|
3119 |
|
|
3120 |
|
|
3121 /* ========================================================================== */ |
|
|
3122 /* === debug_structures ===================================================== */ |
|
|
3123 /* ========================================================================== */ |
|
|
3124 |
|
|
3125 /* |
|
|
3126 At this point, all empty rows and columns are dead. All live columns |
|
|
3127 are "clean" (containing no dead rows) and simplicial (no supercolumns |
|
|
3128 yet). Rows may contain dead columns, but all live rows contain at |
|
|
3129 least one live column. |
|
|
3130 */ |
|
|
3131 |
|
|
3132 PRIVATE void debug_structures |
|
|
3133 ( |
|
|
3134 /* === Parameters ======================================================= */ |
|
|
3135 |
|
|
3136 int n_row, |
|
|
3137 int n_col, |
|
|
3138 Colamd_Row Row [], |
|
|
3139 Colamd_Col Col [], |
|
|
3140 int A [], |
|
|
3141 int n_col2 |
|
|
3142 ) |
|
|
3143 { |
|
|
3144 /* === Local variables ================================================== */ |
|
|
3145 |
|
|
3146 int i ; |
|
|
3147 int c ; |
|
|
3148 int *cp ; |
|
|
3149 int *cp_end ; |
|
|
3150 int len ; |
|
|
3151 int score ; |
|
|
3152 int r ; |
|
|
3153 int *rp ; |
|
|
3154 int *rp_end ; |
|
|
3155 int deg ; |
|
|
3156 |
|
|
3157 /* === Check A, Row, and Col ============================================ */ |
|
|
3158 |
|
|
3159 for (c = 0 ; c < n_col ; c++) |
|
|
3160 { |
|
|
3161 if (COL_IS_ALIVE (c)) |
|
|
3162 { |
|
|
3163 len = Col [c].length ; |
|
|
3164 score = Col [c].shared2.score ; |
|
|
3165 DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; |
|
|
3166 ASSERT (len > 0) ; |
|
|
3167 ASSERT (score >= 0) ; |
|
|
3168 ASSERT (Col [c].shared1.thickness == 1) ; |
|
|
3169 cp = &A [Col [c].start] ; |
|
|
3170 cp_end = cp + len ; |
|
|
3171 while (cp < cp_end) |
|
|
3172 { |
|
|
3173 r = *cp++ ; |
|
|
3174 ASSERT (ROW_IS_ALIVE (r)) ; |
|
|
3175 } |
|
|
3176 } |
|
|
3177 else |
|
|
3178 { |
|
|
3179 i = Col [c].shared2.order ; |
|
|
3180 ASSERT (i >= n_col2 && i < n_col) ; |
|
|
3181 } |
|
|
3182 } |
|
|
3183 |
|
|
3184 for (r = 0 ; r < n_row ; r++) |
|
|
3185 { |
|
|
3186 if (ROW_IS_ALIVE (r)) |
|
|
3187 { |
|
|
3188 i = 0 ; |
|
|
3189 len = Row [r].length ; |
|
|
3190 deg = Row [r].shared1.degree ; |
|
|
3191 ASSERT (len > 0) ; |
|
|
3192 ASSERT (deg > 0) ; |
|
|
3193 rp = &A [Row [r].start] ; |
|
|
3194 rp_end = rp + len ; |
|
|
3195 while (rp < rp_end) |
|
|
3196 { |
|
|
3197 c = *rp++ ; |
|
|
3198 if (COL_IS_ALIVE (c)) |
|
|
3199 { |
|
|
3200 i++ ; |
|
|
3201 } |
|
|
3202 } |
|
|
3203 ASSERT (i > 0) ; |
|
|
3204 } |
|
|
3205 } |
|
|
3206 } |
|
|
3207 |
|
|
3208 |
|
|
3209 /* ========================================================================== */ |
|
|
3210 /* === debug_deg_lists ====================================================== */ |
|
|
3211 /* ========================================================================== */ |
|
|
3212 |
|
|
3213 /* |
|
|
3214 Prints the contents of the degree lists. Counts the number of columns |
|
|
3215 in the degree list and compares it to the total it should have. Also |
|
|
3216 checks the row degrees. |
|
|
3217 */ |
|
|
3218 |
|
|
3219 PRIVATE void debug_deg_lists |
|
|
3220 ( |
|
|
3221 /* === Parameters ======================================================= */ |
|
|
3222 |
|
|
3223 int n_row, |
|
|
3224 int n_col, |
|
|
3225 Colamd_Row Row [], |
|
|
3226 Colamd_Col Col [], |
|
|
3227 int head [], |
|
|
3228 int min_score, |
|
|
3229 int should, |
|
|
3230 int max_deg |
|
|
3231 ) |
|
|
3232 { |
|
|
3233 /* === Local variables ================================================== */ |
|
|
3234 |
|
|
3235 int deg ; |
|
|
3236 int col ; |
|
|
3237 int have ; |
|
|
3238 int row ; |
|
|
3239 |
|
|
3240 /* === Check the degree lists =========================================== */ |
|
|
3241 |
|
|
3242 if (n_col > 10000 && colamd_debug <= 0) |
|
|
3243 { |
|
|
3244 return ; |
|
|
3245 } |
|
|
3246 have = 0 ; |
|
|
3247 DEBUG4 (("Degree lists: %d\n", min_score)) ; |
|
|
3248 for (deg = 0 ; deg <= n_col ; deg++) |
|
|
3249 { |
|
|
3250 col = head [deg] ; |
|
|
3251 if (col == EMPTY) |
|
|
3252 { |
|
|
3253 continue ; |
|
|
3254 } |
|
|
3255 DEBUG4 (("%d:", deg)) ; |
|
|
3256 while (col != EMPTY) |
|
|
3257 { |
|
|
3258 DEBUG4 ((" %d", col)) ; |
|
|
3259 have += Col [col].shared1.thickness ; |
|
|
3260 ASSERT (COL_IS_ALIVE (col)) ; |
|
|
3261 col = Col [col].shared4.degree_next ; |
|
|
3262 } |
|
|
3263 DEBUG4 (("\n")) ; |
|
|
3264 } |
|
|
3265 DEBUG4 (("should %d have %d\n", should, have)) ; |
|
|
3266 ASSERT (should == have) ; |
|
|
3267 |
|
|
3268 /* === Check the row degrees ============================================ */ |
|
|
3269 |
|
|
3270 if (n_row > 10000 && colamd_debug <= 0) |
|
|
3271 { |
|
|
3272 return ; |
|
|
3273 } |
|
|
3274 for (row = 0 ; row < n_row ; row++) |
|
|
3275 { |
|
|
3276 if (ROW_IS_ALIVE (row)) |
|
|
3277 { |
|
|
3278 ASSERT (Row [row].shared1.degree <= max_deg) ; |
|
|
3279 } |
|
|
3280 } |
|
|
3281 } |
|
|
3282 |
|
|
3283 |
|
|
3284 /* ========================================================================== */ |
|
|
3285 /* === debug_mark =========================================================== */ |
|
|
3286 /* ========================================================================== */ |
|
|
3287 |
|
|
3288 /* |
|
|
3289 Ensures that the tag_mark is less that the maximum and also ensures that |
|
|
3290 each entry in the mark array is less than the tag mark. |
|
|
3291 */ |
|
|
3292 |
|
|
3293 PRIVATE void debug_mark |
|
|
3294 ( |
|
|
3295 /* === Parameters ======================================================= */ |
|
|
3296 |
|
|
3297 int n_row, |
|
|
3298 Colamd_Row Row [], |
|
|
3299 int tag_mark, |
|
|
3300 int max_mark |
|
|
3301 ) |
|
|
3302 { |
|
|
3303 /* === Local variables ================================================== */ |
|
|
3304 |
|
|
3305 int r ; |
|
|
3306 |
|
|
3307 /* === Check the Row marks ============================================== */ |
|
|
3308 |
|
|
3309 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; |
|
|
3310 if (n_row > 10000 && colamd_debug <= 0) |
|
|
3311 { |
|
|
3312 return ; |
|
|
3313 } |
|
|
3314 for (r = 0 ; r < n_row ; r++) |
|
|
3315 { |
|
|
3316 ASSERT (Row [r].shared2.mark < tag_mark) ; |
|
|
3317 } |
|
|
3318 } |
|
|
3319 |
|
|
3320 |
|
|
3321 /* ========================================================================== */ |
|
|
3322 /* === debug_matrix ========================================================= */ |
|
|
3323 /* ========================================================================== */ |
|
|
3324 |
|
|
3325 /* |
|
|
3326 Prints out the contents of the columns and the rows. |
|
|
3327 */ |
|
|
3328 |
|
|
3329 PRIVATE void debug_matrix |
|
|
3330 ( |
|
|
3331 /* === Parameters ======================================================= */ |
|
|
3332 |
|
|
3333 int n_row, |
|
|
3334 int n_col, |
|
|
3335 Colamd_Row Row [], |
|
|
3336 Colamd_Col Col [], |
|
|
3337 int A [] |
|
|
3338 ) |
|
|
3339 { |
|
|
3340 /* === Local variables ================================================== */ |
|
|
3341 |
|
|
3342 int r ; |
|
|
3343 int c ; |
|
|
3344 int *rp ; |
|
|
3345 int *rp_end ; |
|
|
3346 int *cp ; |
|
|
3347 int *cp_end ; |
|
|
3348 |
|
|
3349 /* === Dump the rows and columns of the matrix ========================== */ |
|
|
3350 |
|
|
3351 if (colamd_debug < 3) |
|
|
3352 { |
|
|
3353 return ; |
|
|
3354 } |
|
|
3355 DEBUG3 (("DUMP MATRIX:\n")) ; |
|
|
3356 for (r = 0 ; r < n_row ; r++) |
|
|
3357 { |
|
|
3358 DEBUG3 (("Row %d alive? %d\n", r, ROW_IS_ALIVE (r))) ; |
|
|
3359 if (ROW_IS_DEAD (r)) |
|
|
3360 { |
|
|
3361 continue ; |
|
|
3362 } |
|
|
3363 DEBUG3 (("start %d length %d degree %d\n", |
|
|
3364 Row [r].start, Row [r].length, Row [r].shared1.degree)) ; |
|
|
3365 rp = &A [Row [r].start] ; |
|
|
3366 rp_end = rp + Row [r].length ; |
|
|
3367 while (rp < rp_end) |
|
|
3368 { |
|
|
3369 c = *rp++ ; |
|
|
3370 DEBUG4 ((" %d col %d\n", COL_IS_ALIVE (c), c)) ; |
|
|
3371 } |
|
|
3372 } |
|
|
3373 |
|
|
3374 for (c = 0 ; c < n_col ; c++) |
|
|
3375 { |
|
|
3376 DEBUG3 (("Col %d alive? %d\n", c, COL_IS_ALIVE (c))) ; |
|
|
3377 if (COL_IS_DEAD (c)) |
|
|
3378 { |
|
|
3379 continue ; |
|
|
3380 } |
|
|
3381 DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", |
|
|
3382 Col [c].start, Col [c].length, |
|
|
3383 Col [c].shared1.thickness, Col [c].shared2.score)) ; |
|
|
3384 cp = &A [Col [c].start] ; |
|
|
3385 cp_end = cp + Col [c].length ; |
|
|
3386 while (cp < cp_end) |
|
|
3387 { |
|
|
3388 r = *cp++ ; |
|
|
3389 DEBUG4 ((" %d row %d\n", ROW_IS_ALIVE (r), r)) ; |
|
|
3390 } |
|
|
3391 } |
|
|
3392 } |
|
|
3393 |
|
|
3394 PRIVATE void colamd_get_debug |
|
|
3395 ( |
|
|
3396 char *method |
|
|
3397 ) |
|
|
3398 { |
|
|
3399 colamd_debug = 0 ; /* no debug printing */ |
|
|
3400 |
|
|
3401 /* get "D" environment variable, which gives the debug printing level */ |
|
|
3402 if (getenv ("D")) |
|
|
3403 { |
|
|
3404 colamd_debug = atoi (getenv ("D")) ; |
|
|
3405 } |
|
|
3406 |
|
|
3407 DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", |
|
|
3408 method, colamd_debug)) ; |
|
|
3409 } |
|
|
3410 |
|
|
3411 #endif /* NDEBUG */ |
|
|
3412 |
