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