1 /* ========================================================================== */
2 /* === CCOLAMD/CSYMAMD - a constrained column ordering algorithm ============ */
3 /* ========================================================================== */
4
5 /* ----------------------------------------------------------------------------
6 * CCOLAMD, Copyright (C) Univ. of Florida. Authors: Timothy A. Davis,
7 * Sivasankaran Rajamanickam, and Stefan Larimore
8 * -------------------------------------------------------------------------- */
9
10 /*
11 * ccolamd: a constrained approximate minimum degree column ordering
12 * algorithm, LU factorization of symmetric or unsymmetric matrices,
13 * QR factorization, least squares, interior point methods for
14 * linear programming problems, and other related problems.
15 *
16 * csymamd: a constrained approximate minimum degree ordering algorithm for
17 * Cholesky factorization of symmetric matrices.
18 *
19 * Purpose:
20 *
21 * CCOLAMD computes a permutation Q such that the Cholesky factorization of
22 * (AQ)'(AQ) has less fill-in and requires fewer floating point operations
23 * than A'A. This also provides a good ordering for sparse partial
24 * pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
25 * factorization, and P is computed during numerical factorization via
26 * conventional partial pivoting with row interchanges. CCOLAMD is an
27 * extension of COLAMD, available as built-in function in MATLAB Version 6,
28 * available from MathWorks, Inc. (http://www.mathworks.com). This
29 * routine can be used in place of COLAMD in MATLAB.
30 *
31 * CSYMAMD computes a permutation P of a symmetric matrix A such that the
32 * Cholesky factorization of PAP' has less fill-in and requires fewer
33 * floating point operations than A. CSYMAMD constructs a matrix M such
34 * that M'M has the same nonzero pattern of A, and then orders the columns
35 * of M using colmmd. The column ordering of M is then returned as the
36 * row and column ordering P of A. CSYMAMD is an extension of SYMAMD.
37 *
38 * Authors:
39 *
40 * Timothy A. Davis and S. Rajamanickam wrote CCOLAMD, based directly on
41 * COLAMD by Stefan I. Larimore and Timothy A. Davis, University of
42 * Florida. The algorithm was developed in collaboration with John
43 * Gilbert, (UCSB, then at Xerox PARC), and Esmond Ng, (Lawrence Berkeley
44 * National Lab, then at Oak Ridge National Laboratory).
45 *
46 * Acknowledgements:
47 *
48 * This work was supported by the National Science Foundation, under
49 * grants DMS-9504974 and DMS-9803599, CCR-0203270, and a grant from the
50 * Sandia National Laboratory (Dept. of Energy).
51 *
52 * Copyright and License:
53 *
54 * Copyright (c) 1998-2005 by the University of Florida.
55 * All Rights Reserved.
56 * COLAMD is also available under alternate licenses, contact T. Davis
57 * for details.
58 *
59 * See CCOLAMD/Doc/License.txt for the license.
60 *
61 * Availability:
62 *
63 * The CCOLAMD/CSYMAMD library is available at
64 *
65 * http://www.suitesparse.com
66 *
67 * See the ChangeLog file for changes since Version 1.0.
68 */
69
70 /* ========================================================================== */
71 /* === Description of user-callable routines ================================ */
72 /* ========================================================================== */
73
74 /* CCOLAMD includes both int and SuiteSparse_long versions of all its routines.
75 * The description below is for the int version. For SuiteSparse_long, all
76 * int arguments become SuiteSparse_long integers. SuiteSparse_long is
77 * normally defined as long, except for WIN64 */
78
79 /* ----------------------------------------------------------------------------
80 * ccolamd_recommended:
81 * ----------------------------------------------------------------------------
82 *
83 * C syntax:
84 *
85 * #include "ccolamd.h"
86 * size_t ccolamd_recommended (int nnz, int n_row, int n_col) ;
87 * size_t ccolamd_l_recommended (SuiteSparse_long nnz,
88 * SuiteSparse_long n_row, SuiteSparse_long n_col) ;
89 *
90 * Purpose:
91 *
92 * Returns recommended value of Alen for use by ccolamd. Returns 0
93 * if any input argument is negative. The use of this routine
94 * is optional. Not needed for csymamd, which dynamically allocates
95 * its own memory.
96 *
97 * Arguments (all input arguments):
98 *
99 * int nnz ; Number of nonzeros in the matrix A. This must
100 * be the same value as p [n_col] in the call to
101 * ccolamd - otherwise you will get a wrong value
102 * of the recommended memory to use.
103 *
104 * int n_row ; Number of rows in the matrix A.
105 *
106 * int n_col ; Number of columns in the matrix A.
107 *
108 * ----------------------------------------------------------------------------
109 * ccolamd_set_defaults:
110 * ----------------------------------------------------------------------------
111 *
112 * C syntax:
113 *
114 * #include "ccolamd.h"
115 * ccolamd_set_defaults (double knobs [CCOLAMD_KNOBS]) ;
116 * ccolamd_l_set_defaults (double knobs [CCOLAMD_KNOBS]) ;
117 *
118 * Purpose:
119 *
120 * Sets the default parameters. The use of this routine is optional.
121 * Passing a (double *) NULL pointer for the knobs results in the
122 * default parameter settings.
123 *
124 * Arguments:
125 *
126 * double knobs [CCOLAMD_KNOBS] ; Output only.
127 *
128 * knobs [0] and knobs [1] behave differently than they did in COLAMD.
129 * The other knobs are new to CCOLAMD.
130 *
131 * knobs [0]: dense row control
132 *
133 * For CCOLAMD, rows with more than
134 * max (16, knobs [CCOLAMD_DENSE_ROW] * sqrt (n_col))
135 * entries are removed prior to ordering.
136 *
137 * For CSYMAMD, rows and columns with more than
138 * max (16, knobs [CCOLAMD_DENSE_ROW] * sqrt (n))
139 * entries are removed prior to ordering, and placed last in the
140 * output ordering (subject to the constraints).
141 *
142 * If negative, only completely dense rows are removed. If you
143 * intend to use CCOLAMD for a Cholesky factorization of A*A', set
144 * knobs [CCOLAMD_DENSE_ROW] to -1, which is more appropriate for
145 * that case.
146 *
147 * Default: 10.
148 *
149 * knobs [1]: dense column control
150 *
151 * For CCOLAMD, columns with more than
152 * max (16, knobs [CCOLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
153 * entries are removed prior to ordering, and placed last in the
154 * output column ordering (subject to the constraints).
155 * Not used by CSYMAMD. If negative, only completely dense
156 * columns are removed. Default: 10.
157 *
158 * knobs [2]: aggressive absorption
159 *
160 * knobs [CCOLAMD_AGGRESSIVE] controls whether or not to do
161 * aggressive absorption during the ordering. Default is TRUE
162 * (nonzero). If zero, no aggressive absorption is performed.
163 *
164 * knobs [3]: optimize ordering for LU or Cholesky
165 *
166 * knobs [CCOLAMD_LU] controls an option that optimizes the
167 * ordering for the LU of A or the Cholesky factorization of A'A.
168 * If TRUE (nonzero), an ordering optimized for LU is performed.
169 * If FALSE (zero), an ordering for Cholesky is performed.
170 * Default is FALSE. CSYMAMD ignores this parameter; it always
171 * orders for Cholesky.
172 *
173 * ----------------------------------------------------------------------------
174 * ccolamd:
175 * ----------------------------------------------------------------------------
176 *
177 * C syntax:
178 *
179 * #include "ccolamd.h"
180 * int ccolamd (int n_row, int n_col, int Alen, int *A, int *p,
181 * double knobs [CCOLAMD_KNOBS], int stats [CCOLAMD_STATS],
182 * int *cmember) ;
183 *
184 * SuiteSparse_long ccolamd_l (SuiteSparse_long n_row,
185 * SuiteSparse_long n_col, SuiteSparse_long Alen,
186 * SuiteSparse_long *A, SuiteSparse_long *p,
187 * double knobs [CCOLAMD_KNOBS],
188 * SuiteSparse_long stats [CCOLAMD_STATS],
189 * SuiteSparse_long *cmember) ;
190 *
191 * Purpose:
192 *
193 * Computes a column ordering (Q) of A such that P(AQ)=LU or
194 * (AQ)'AQ=LL' have less fill-in and require fewer floating point
195 * operations than factorizing the unpermuted matrix A or A'A,
196 * respectively.
197 *
198 * Returns:
199 *
200 * TRUE (1) if successful, FALSE (0) otherwise.
201 *
202 * Arguments (for int version):
203 *
204 * int n_row ; Input argument.
205 *
206 * Number of rows in the matrix A.
207 * Restriction: n_row >= 0.
208 * ccolamd returns FALSE if n_row is negative.
209 *
210 * int n_col ; Input argument.
211 *
212 * Number of columns in the matrix A.
213 * Restriction: n_col >= 0.
214 * ccolamd returns FALSE if n_col is negative.
215 *
216 * int Alen ; Input argument.
217 *
218 * Restriction (see note):
219 * Alen >= MAX (2*nnz, 4*n_col) + 17*n_col + 7*n_row + 7, where
220 * nnz = p [n_col]. ccolamd returns FALSE if this condition is
221 * not met. We recommend about nnz/5 more space for better
222 * efficiency. This restriction makes an modest assumption
223 * regarding the size of two typedef'd structures in ccolamd.h.
224 * We do, however, guarantee that
225 *
226 * Alen >= ccolamd_recommended (nnz, n_row, n_col)
227 *
228 * will work efficiently.
229 *
230 * int A [Alen] ; Input argument, undefined on output.
231 *
232 * A is an integer array of size Alen. Alen must be at least as
233 * large as the bare minimum value given above, but this is very
234 * low, and can result in excessive run time. For best
235 * performance, we recommend that Alen be greater than or equal to
236 * ccolamd_recommended (nnz, n_row, n_col), which adds
237 * nnz/5 to the bare minimum value given above.
238 *
239 * On input, the row indices of the entries in column c of the
240 * matrix are held in A [(p [c]) ... (p [c+1]-1)]. The row indices
241 * in a given column c need not be in ascending order, and
242 * duplicate row indices may be be present. However, ccolamd will
243 * work a little faster if both of these conditions are met
244 * (ccolamd puts the matrix into this format, if it finds that the
245 * the conditions are not met).
246 *
247 * The matrix is 0-based. That is, rows are in the range 0 to
248 * n_row-1, and columns are in the range 0 to n_col-1. ccolamd
249 * returns FALSE if any row index is out of range.
250 *
251 * The contents of A are modified during ordering, and are
252 * undefined on output.
253 *
254 * int p [n_col+1] ; Both input and output argument.
255 *
256 * p is an integer array of size n_col+1. On input, it holds the
257 * "pointers" for the column form of the matrix A. Column c of
258 * the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
259 * entry, p [0], must be zero, and p [c] <= p [c+1] must hold
260 * for all c in the range 0 to n_col-1. The value nnz = p [n_col]
261 * is thus the total number of entries in the pattern of the
262 * matrix A. ccolamd returns FALSE if these conditions are not
263 * met.
264 *
265 * On output, if ccolamd returns TRUE, the array p holds the column
266 * permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
267 * the first column index in the new ordering, and p [n_col-1] is
268 * the last. That is, p [k] = j means that column j of A is the
269 * kth pivot column, in AQ, where k is in the range 0 to n_col-1
270 * (p [0] = j means that column j of A is the first column in AQ).
271 *
272 * If ccolamd returns FALSE, then no permutation is returned, and
273 * p is undefined on output.
274 *
275 * double knobs [CCOLAMD_KNOBS] ; Input argument.
276 *
277 * See ccolamd_set_defaults for a description.
278 *
279 * int stats [CCOLAMD_STATS] ; Output argument.
280 *
281 * Statistics on the ordering, and error status.
282 * See ccolamd.h for related definitions.
283 * ccolamd returns FALSE if stats is not present.
284 *
285 * stats [0]: number of dense or empty rows ignored.
286 *
287 * stats [1]: number of dense or empty columns ignored (and
288 * ordered last in the output permutation p, subject to the
289 * constraints). Note that a row can become "empty" if it
290 * contains only "dense" and/or "empty" columns, and similarly
291 * a column can become "empty" if it only contains "dense"
292 * and/or "empty" rows.
293 *
294 * stats [2]: number of garbage collections performed. This can
295 * be excessively high if Alen is close to the minimum
296 * required value.
297 *
298 * stats [3]: status code. < 0 is an error code.
299 * > 1 is a warning or notice.
300 *
301 * 0 OK. Each column of the input matrix contained row
302 * indices in increasing order, with no duplicates.
303 *
304 * 1 OK, but columns of input matrix were jumbled (unsorted
305 * columns or duplicate entries). CCOLAMD had to do some
306 * extra work to sort the matrix first and remove
307 * duplicate entries, but it still was able to return a
308 * valid permutation (return value of ccolamd was TRUE).
309 *
310 * stats [4]: highest column index of jumbled columns
311 * stats [5]: last seen duplicate or unsorted row index
312 * stats [6]: number of duplicate or unsorted row indices
313 *
314 * -1 A is a null pointer
315 *
316 * -2 p is a null pointer
317 *
318 * -3 n_row is negative. stats [4]: n_row
319 *
320 * -4 n_col is negative. stats [4]: n_col
321 *
322 * -5 number of nonzeros in matrix is negative
323 *
324 * stats [4]: number of nonzeros, p [n_col]
325 *
326 * -6 p [0] is nonzero
327 *
328 * stats [4]: p [0]
329 *
330 * -7 A is too small
331 *
332 * stats [4]: required size
333 * stats [5]: actual size (Alen)
334 *
335 * -8 a column has a negative number of entries
336 *
337 * stats [4]: column with < 0 entries
338 * stats [5]: number of entries in col
339 *
340 * -9 a row index is out of bounds
341 *
342 * stats [4]: column with bad row index
343 * stats [5]: bad row index
344 * stats [6]: n_row, # of rows of matrx
345 *
346 * -10 (unused; see csymamd)
347 *
348 * int cmember [n_col] ; Input argument.
349 *
350 * cmember is new to CCOLAMD. It did not appear in COLAMD.
351 * It places contraints on the output ordering. s = cmember [j]
352 * gives the constraint set s that contains the column j
353 * (Restriction: 0 <= s < n_col). In the output column
354 * permutation, all columns in set 0 appear first, followed by
355 * all columns in set 1, and so on. If NULL, all columns are
356 * treated as if they were in a single constraint set, and you
357 * will obtain the same ordering as COLAMD (with one exception:
358 * the dense row/column threshold and other default knobs in
359 * CCOLAMD and COLAMD are different).
360 *
361 * Example:
362 *
363 * See ccolamd_example.c for a complete example.
364 *
365 * To order the columns of a 5-by-4 matrix with 11 nonzero entries in
366 * the following nonzero pattern
367 *
368 * x 0 x 0
369 * x 0 x x
370 * 0 x x 0
371 * 0 0 x x
372 * x x 0 0
373 *
374 * with default knobs, no output statistics, and no ordering
375 * constraints, do the following:
376 *
377 * #include "ccolamd.h"
378 * #define ALEN 144
379 * int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
380 * int p [ ] = {0, 3, 5, 9, 11} ;
381 * int stats [CCOLAMD_STATS] ;
382 * ccolamd (5, 4, ALEN, A, p, (double *) NULL, stats, NULL) ;
383 *
384 * The permutation is returned in the array p, and A is destroyed.
385 *
386 * ----------------------------------------------------------------------------
387 * csymamd:
388 * ----------------------------------------------------------------------------
389 *
390 * C syntax:
391 *
392 * #include "ccolamd.h"
393 *
394 * int csymamd (int n, int *A, int *p, int *perm,
395 * double knobs [CCOLAMD_KNOBS], int stats [CCOLAMD_STATS],
396 * void (*allocate) (size_t, size_t), void (*release) (void *),
397 * int *cmember, int stype) ;
398 *
399 * SuiteSparse_long csymamd_l (SuiteSparse_long n,
400 * SuiteSparse_long *A, SuiteSparse_long *p,
401 * SuiteSparse_long *perm, double knobs [CCOLAMD_KNOBS],
402 * SuiteSparse_long stats [CCOLAMD_STATS], void (*allocate)
403 * (size_t, size_t), void (*release) (void *),
404 * SuiteSparse_long *cmember, SuiteSparse_long stype) ;
405 *
406 * Purpose:
407 *
408 * The csymamd routine computes an ordering P of a symmetric sparse
409 * matrix A such that the Cholesky factorization PAP' = LL' remains
410 * sparse. It is based on a column ordering of a matrix M constructed
411 * so that the nonzero pattern of M'M is the same as A. Either the
412 * lower or upper triangular part of A can be used, or the pattern
413 * A+A' can be used. You must pass your selected memory allocator
414 * (usually calloc/free or mxCalloc/mxFree) to csymamd, for it to
415 * allocate memory for the temporary matrix M.
416 *
417 * Returns:
418 *
419 * TRUE (1) if successful, FALSE (0) otherwise.
420 *
421 * Arguments:
422 *
423 * int n ; Input argument.
424 *
425 * Number of rows and columns in the symmetrix matrix A.
426 * Restriction: n >= 0.
427 * csymamd returns FALSE if n is negative.
428 *
429 * int A [nnz] ; Input argument.
430 *
431 * A is an integer array of size nnz, where nnz = p [n].
432 *
433 * The row indices of the entries in column c of the matrix are
434 * held in A [(p [c]) ... (p [c+1]-1)]. The row indices in a
435 * given column c need not be in ascending order, and duplicate
436 * row indices may be present. However, csymamd will run faster
437 * if the columns are in sorted order with no duplicate entries.
438 *
439 * The matrix is 0-based. That is, rows are in the range 0 to
440 * n-1, and columns are in the range 0 to n-1. csymamd
441 * returns FALSE if any row index is out of range.
442 *
443 * The contents of A are not modified.
444 *
445 * int p [n+1] ; Input argument.
446 *
447 * p is an integer array of size n+1. On input, it holds the
448 * "pointers" for the column form of the matrix A. Column c of
449 * the matrix A is held in A [(p [c]) ... (p [c+1]-1)]. The first
450 * entry, p [0], must be zero, and p [c] <= p [c+1] must hold
451 * for all c in the range 0 to n-1. The value p [n] is
452 * thus the total number of entries in the pattern of the matrix A.
453 * csymamd returns FALSE if these conditions are not met.
454 *
455 * The contents of p are not modified.
456 *
457 * int perm [n+1] ; Output argument.
458 *
459 * On output, if csymamd returns TRUE, the array perm holds the
460 * permutation P, where perm [0] is the first index in the new
461 * ordering, and perm [n-1] is the last. That is, perm [k] = j
462 * means that row and column j of A is the kth column in PAP',
463 * where k is in the range 0 to n-1 (perm [0] = j means
464 * that row and column j of A are the first row and column in
465 * PAP'). The array is used as a workspace during the ordering,
466 * which is why it must be of length n+1, not just n.
467 *
468 * double knobs [CCOLAMD_KNOBS] ; Input argument.
469 *
470 * See colamd_set_defaults for a description.
471 *
472 * int stats [CCOLAMD_STATS] ; Output argument.
473 *
474 * Statistics on the ordering, and error status.
475 * See ccolamd.h for related definitions.
476 * csymand returns FALSE if stats is not present.
477 *
478 * stats [0]: number of dense or empty row and columns ignored
479 * (and ordered last in the output permutation perm, subject
480 * to the constraints). Note that a row/column can become
481 * "empty" if it contains only "dense" and/or "empty"
482 * columns/rows.
483 *
484 * stats [1]: (same as stats [0])
485 *
486 * stats [2]: number of garbage collections performed.
487 *
488 * stats [3]: status code. < 0 is an error code.
489 * > 1 is a warning or notice.
490 *
491 * 0 to -9: same as ccolamd, with n replacing n_col and n_row,
492 * and -3 and -7 are unused.
493 *
494 * -10 out of memory (unable to allocate temporary workspace
495 * for M or count arrays using the "allocate" routine
496 * passed into csymamd).
497 *
498 * void * (*allocate) (size_t, size_t)
499 *
500 * A pointer to a function providing memory allocation. The
501 * allocated memory must be returned initialized to zero. For a
502 * C application, this argument should normally be a pointer to
503 * calloc. For a MATLAB mexFunction, the routine mxCalloc is
504 * passed instead.
505 *
506 * void (*release) (size_t, size_t)
507 *
508 * A pointer to a function that frees memory allocated by the
509 * memory allocation routine above. For a C application, this
510 * argument should normally be a pointer to free. For a MATLAB
511 * mexFunction, the routine mxFree is passed instead.
512 *
513 * int cmember [n] ; Input argument.
514 *
515 * Same as ccolamd, except that cmember is of size n, and it places
516 * contraints symmetrically, on both the row and column ordering.
517 * Entries in cmember must be in the range 0 to n-1.
518 *
519 * int stype ; Input argument.
520 *
521 * If stype < 0, then only the strictly lower triangular part of
522 * A is accessed. The upper triangular part is assumed to be the
523 * transpose of the lower triangular part. This is the same as
524 * SYMAMD, which did not have an stype parameter.
525 *
526 * If stype > 0, only the strictly upper triangular part of A is
527 * accessed. The lower triangular part is assumed to be the
528 * transpose of the upper triangular part.
529 *
530 * If stype == 0, then the nonzero pattern of A+A' is ordered.
531 *
532 * ----------------------------------------------------------------------------
533 * ccolamd_report:
534 * ----------------------------------------------------------------------------
535 *
536 * C syntax:
537 *
538 * #include "ccolamd.h"
539 * ccolamd_report (int stats [CCOLAMD_STATS]) ;
540 * ccolamd_l_report (SuiteSparse_long stats [CCOLAMD_STATS]) ;
541 *
542 * Purpose:
543 *
544 * Prints the error status and statistics recorded in the stats
545 * array on the standard error output (for a standard C routine)
546 * or on the MATLAB output (for a mexFunction).
547 *
548 * Arguments:
549 *
550 * int stats [CCOLAMD_STATS] ; Input only. Statistics from ccolamd.
551 *
552 *
553 * ----------------------------------------------------------------------------
554 * csymamd_report:
555 * ----------------------------------------------------------------------------
556 *
557 * C syntax:
558 *
559 * #include "ccolamd.h"
560 * csymamd_report (int stats [CCOLAMD_STATS]) ;
561 * csymamd_l_report (SuiteSparse_long stats [CCOLAMD_STATS]) ;
562 *
563 * Purpose:
564 *
565 * Prints the error status and statistics recorded in the stats
566 * array on the standard error output (for a standard C routine)
567 * or on the MATLAB output (for a mexFunction).
568 *
569 * Arguments:
570 *
571 * int stats [CCOLAMD_STATS] ; Input only. Statistics from csymamd.
572 *
573 */
574
575
576 /* ========================================================================== */
577 /* === Scaffolding code definitions ======================================== */
578 /* ========================================================================== */
579
580 /* Ensure that debugging is turned off: */
581 #ifndef NDEBUG
582 #define NDEBUG
583 #endif
584
585 /* turn on debugging by uncommenting the following line
586 #undef NDEBUG
587 */
588
589 /* ========================================================================== */
590 /* === Include files ======================================================== */
591 /* ========================================================================== */
592
593 #include "ccolamd.h"
594
595 #include <stdlib.h>
596 #include <math.h>
597 #include <limits.h>
598
599 #ifdef MATLAB_MEX_FILE
600 #include "mex.h"
601 #include "matrix.h"
602 #endif
603
604 #if !defined (NPRINT) || !defined (NDEBUG)
605 #include <stdio.h>
606 #endif
607
608 #ifndef NULL
609 #define NULL ((void *) 0)
610 #endif
611
612 /* ========================================================================== */
613 /* === int or SuiteSparse_long ============================================== */
614 /* ========================================================================== */
615
616 #ifdef DLONG
617
618 #define Int SuiteSparse_long
619 #define ID SuiteSparse_long_id
620 #define Int_MAX SuiteSparse_long_max
621
622 #define CCOLAMD_recommended ccolamd_l_recommended
623 #define CCOLAMD_set_defaults ccolamd_l_set_defaults
624 #define CCOLAMD_2 ccolamd2_l
625 #define CCOLAMD_MAIN ccolamd_l
626 #define CCOLAMD_apply_order ccolamd_l_apply_order
627 #define CCOLAMD_postorder ccolamd_l_postorder
628 #define CCOLAMD_post_tree ccolamd_l_post_tree
629 #define CCOLAMD_fsize ccolamd_l_fsize
630 #define CSYMAMD_MAIN csymamd_l
631 #define CCOLAMD_report ccolamd_l_report
632 #define CSYMAMD_report csymamd_l_report
633
634 #else
635
636 #define Int int
637 #define ID "%d"
638 #define Int_MAX INT_MAX
639
640 #define CCOLAMD_recommended ccolamd_recommended
641 #define CCOLAMD_set_defaults ccolamd_set_defaults
642 #define CCOLAMD_2 ccolamd2
643 #define CCOLAMD_MAIN ccolamd
644 #define CCOLAMD_apply_order ccolamd_apply_order
645 #define CCOLAMD_postorder ccolamd_postorder
646 #define CCOLAMD_post_tree ccolamd_post_tree
647 #define CCOLAMD_fsize ccolamd_fsize
648 #define CSYMAMD_MAIN csymamd
649 #define CCOLAMD_report ccolamd_report
650 #define CSYMAMD_report csymamd_report
651
652 #endif
653
654 /* ========================================================================== */
655 /* === Row and Column structures ============================================ */
656 /* ========================================================================== */
657
658 typedef struct CColamd_Col_struct
659 {
660 /* size of this struct is 8 integers if no padding occurs */
661
662 Int start ; /* index for A of first row in this column, or DEAD */
663 /* if column is dead */
664 Int length ; /* number of rows in this column */
665 union
666 {
667 Int thickness ; /* number of original columns represented by this */
668 /* col, if the column is alive */
669 Int parent ; /* parent in parent tree super-column structure, if */
670 /* the column is dead */
671 } shared1 ;
672 union
673 {
674 Int score ;
675 Int order ;
676 } shared2 ;
677 union
678 {
679 Int headhash ; /* head of a hash bucket, if col is at the head of */
680 /* a degree list */
681 Int hash ; /* hash value, if col is not in a degree list */
682 Int prev ; /* previous column in degree list, if col is in a */
683 /* degree list (but not at the head of a degree list) */
684 } shared3 ;
685 union
686 {
687 Int degree_next ; /* next column, if col is in a degree list */
688 Int hash_next ; /* next column, if col is in a hash list */
689 } shared4 ;
690
691 Int nextcol ; /* next column in this supercolumn */
692 Int lastcol ; /* last column in this supercolumn */
693
694 } CColamd_Col ;
695
696
697 typedef struct CColamd_Row_struct
698 {
699 /* size of this struct is 6 integers if no padding occurs */
700
701 Int start ; /* index for A of first col in this row */
702 Int length ; /* number of principal columns in this row */
703 union
704 {
705 Int degree ; /* number of principal & non-principal columns in row */
706 Int p ; /* used as a row pointer in init_rows_cols () */
707 } shared1 ;
708 union
709 {
710 Int mark ; /* for computing set differences and marking dead rows*/
711 Int first_column ;/* first column in row (used in garbage collection) */
712 } shared2 ;
713
714 Int thickness ; /* number of original rows represented by this row */
715 /* that are not yet pivotal */
716 Int front ; /* -1 if an original row */
717 /* k if this row represents the kth frontal matrix */
718 /* where k goes from 0 to at most n_col-1 */
719
720 } CColamd_Row ;
721
722 /* ========================================================================== */
723 /* === basic definitions ==================================================== */
724 /* ========================================================================== */
725
726 #define EMPTY (-1)
727 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
728 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
729
730 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
731 #define GLOBAL
732 #define PUBLIC
733 #define PRIVATE static
734
735 #define DENSE_DEGREE(alpha,n) \
736 ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
737
738 #define CMEMBER(c) ((cmember == (Int *) NULL) ? (0) : (cmember [c]))
739
740 /* True if x is NaN */
741 #define SCALAR_IS_NAN(x) ((x) != (x))
742
743 /* true if an integer (stored in double x) would overflow (or if x is NaN) */
744 #define INT_OVERFLOW(x) ((!((x) * (1.0+1e-8) <= (double) Int_MAX)) \
745 || SCALAR_IS_NAN (x))
746
747 #define ONES_COMPLEMENT(r) (-(r)-1)
748 #undef TRUE
749 #undef FALSE
750 #define TRUE (1)
751 #define FALSE (0)
752
753 /* Row and column status */
754 #define ALIVE (0)
755 #define DEAD (-1)
756
757 /* Column status */
758 #define DEAD_PRINCIPAL (-1)
759 #define DEAD_NON_PRINCIPAL (-2)
760
761 /* Macros for row and column status update and checking. */
762 #define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
763 #define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE)
764 #define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE)
765 #define COL_IS_DEAD(c) (Col [c].start < ALIVE)
766 #define COL_IS_ALIVE(c) (Col [c].start >= ALIVE)
767 #define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL)
768 #define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; }
769 #define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; }
770 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
771
772
773 /* ========================================================================== */
774 /* === ccolamd reporting mechanism ========================================== */
775 /* ========================================================================== */
776
777 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
778 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
779 #define INDEX(i) ((i)+1)
780 #else
781 /* In C, matrices are 0-based and indices are reported as such in *_report */
782 #define INDEX(i) (i)
783 #endif
784
785
786 /* ========================================================================== */
787 /* === Debugging prototypes and definitions ================================= */
788 /* ========================================================================== */
789
790 #ifndef NDEBUG
791
792 #include <assert.h>
793
794 /* debug print level, present only when debugging */
795 PRIVATE Int ccolamd_debug ;
796
797 /* debug print statements */
798 #define DEBUG0(params) { SUITESPARSE_PRINTF (params) ; }
799 #define DEBUG1(params) { if (ccolamd_debug >= 1) SUITESPARSE_PRINTF (params) ; }
800 #define DEBUG2(params) { if (ccolamd_debug >= 2) SUITESPARSE_PRINTF (params) ; }
801 #define DEBUG3(params) { if (ccolamd_debug >= 3) SUITESPARSE_PRINTF (params) ; }
802 #define DEBUG4(params) { if (ccolamd_debug >= 4) SUITESPARSE_PRINTF (params) ; }
803
804 #ifdef MATLAB_MEX_FILE
805 #define ASSERT(expression) (mxAssert ((expression), ""))
806 #else
807 #define ASSERT(expression) (assert (expression))
808 #endif
809
810 PRIVATE void ccolamd_get_debug
811 (
812 char *method
813 ) ;
814
815 PRIVATE void debug_mark
816 (
817 Int n_row,
818 CColamd_Row Row [],
819 Int tag_mark,
820 Int max_mark
821 ) ;
822
823 PRIVATE void debug_matrix
824 (
825 Int n_row,
826 Int n_col,
827 CColamd_Row Row [],
828 CColamd_Col Col [],
829 Int A []
830 ) ;
831
832 PRIVATE void debug_structures
833 (
834 Int n_row,
835 Int n_col,
836 CColamd_Row Row [],
837 CColamd_Col Col [],
838 Int A [],
839 Int in_cset [],
840 Int cset_start []
841 ) ;
842
843 PRIVATE void dump_super
844 (
845 Int super_c,
846 CColamd_Col Col [],
847 Int n_col
848 ) ;
849
850 PRIVATE void debug_deg_lists
851 (
852 Int n_row,
853 Int n_col,
854 CColamd_Row Row [ ],
855 CColamd_Col Col [ ],
856 Int head [ ],
857 Int min_score,
858 Int should,
859 Int max_deg
860 ) ;
861
862 #else
863
864 /* === No debugging ========================================================= */
865
866 #define DEBUG0(params) ;
867 #define DEBUG1(params) ;
868 #define DEBUG2(params) ;
869 #define DEBUG3(params) ;
870 #define DEBUG4(params) ;
871
872 #define ASSERT(expression)
873
874 #endif
875
876 /* ========================================================================== */
877 /* === Prototypes of PRIVATE routines ======================================= */
878 /* ========================================================================== */
879
880 PRIVATE Int init_rows_cols
881 (
882 Int n_row,
883 Int n_col,
884 CColamd_Row Row [ ],
885 CColamd_Col Col [ ],
886 Int A [ ],
887 Int p [ ],
888 Int stats [CCOLAMD_STATS]
889 ) ;
890
891 PRIVATE void init_scoring
892 (
893 Int n_row,
894 Int n_col,
895 CColamd_Row Row [ ],
896 CColamd_Col Col [ ],
897 Int A [ ],
898 Int head [ ],
899 double knobs [CCOLAMD_KNOBS],
900 Int *p_n_row2,
901 Int *p_n_col2,
902 Int *p_max_deg,
903 Int cmember [ ],
904 Int n_cset,
905 Int cset_start [ ],
906 Int dead_cols [ ],
907 Int *p_ndense_row, /* number of dense rows */
908 Int *p_nempty_row, /* number of original empty rows */
909 Int *p_nnewlyempty_row, /* number of newly empty rows */
910 Int *p_ndense_col, /* number of dense cols (excl "empty" cols) */
911 Int *p_nempty_col, /* number of original empty cols */
912 Int *p_nnewlyempty_col /* number of newly empty cols */
913 ) ;
914
915 PRIVATE Int find_ordering
916 (
917 Int n_row,
918 Int n_col,
919 Int Alen,
920 CColamd_Row Row [ ],
921 CColamd_Col Col [ ],
922 Int A [ ],
923 Int head [ ],
924 #ifndef NDEBUG
925 Int n_col2,
926 #endif
927 Int max_deg,
928 Int pfree,
929 Int cset [ ],
930 Int cset_start [ ],
931 #ifndef NDEBUG
932 Int n_cset,
933 #endif
934 Int cmember [ ],
935 Int Front_npivcol [ ],
936 Int Front_nrows [ ],
937 Int Front_ncols [ ],
938 Int Front_parent [ ],
939 Int Front_cols [ ],
940 Int *p_nfr,
941 Int aggressive,
942 Int InFront [ ],
943 Int order_for_lu
944 ) ;
945
946 PRIVATE void detect_super_cols
947 (
948 #ifndef NDEBUG
949 Int n_col,
950 CColamd_Row Row [ ],
951 #endif
952 CColamd_Col Col [ ],
953 Int A [ ],
954 Int head [ ],
955 Int row_start,
956 Int row_length,
957 Int in_set [ ]
958 ) ;
959
960 PRIVATE Int garbage_collection
961 (
962 Int n_row,
963 Int n_col,
964 CColamd_Row Row [ ],
965 CColamd_Col Col [ ],
966 Int A [ ],
967 Int *pfree
968 ) ;
969
970 PRIVATE Int clear_mark
971 (
972 Int tag_mark,
973 Int max_mark,
974 Int n_row,
975 CColamd_Row Row [ ]
976 ) ;
977
978 PRIVATE void print_report
979 (
980 char *method,
981 Int stats [CCOLAMD_STATS]
982 ) ;
983
984
985 /* ========================================================================== */
986 /* === USER-CALLABLE ROUTINES: ============================================== */
987 /* ========================================================================== */
988
989
990 /* ========================================================================== */
991 /* === ccolamd_recommended ================================================== */
992 /* ========================================================================== */
993
994 /*
995 * The ccolamd_recommended routine returns the suggested size for Alen. This
996 * value has been determined to provide good balance between the number of
997 * garbage collections and the memory requirements for ccolamd. If any
998 * argument is negative, or if integer overflow occurs, a 0 is returned as
999 * an error condition.
1000 *
1001 * 2*nnz space is required for the row and column indices of the matrix
1002 * (or 4*n_col, which ever is larger).
1003 *
1004 * CCOLAMD_C (n_col) + CCOLAMD_R (n_row) space is required for the Col and Row
1005 * arrays, respectively, which are internal to ccolamd. This is equal to
1006 * 8*n_col + 6*n_row if the structures are not padded.
1007 *
1008 * An additional n_col space is the minimal amount of "elbow room",
1009 * and nnz/5 more space is recommended for run time efficiency.
1010 *
1011 * The remaining (((3 * n_col) + 1) + 5 * (n_col + 1) + n_row) space is
1012 * for other workspace used in ccolamd which did not appear in colamd.
1013 */
1014
1015 /* add two values of type size_t, and check for integer overflow */
t_add(size_t a,size_t b,int * ok)1016 static size_t t_add (size_t a, size_t b, int *ok)
1017 {
1018 (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
1019 return ((*ok) ? (a + b) : 0) ;
1020 }
1021
1022 /* compute a*k where k is a small integer, and check for integer overflow */
t_mult(size_t a,size_t k,int * ok)1023 static size_t t_mult (size_t a, size_t k, int *ok)
1024 {
1025 size_t i, s = 0 ;
1026 for (i = 0 ; i < k ; i++)
1027 {
1028 s = t_add (s, a, ok) ;
1029 }
1030 return (s) ;
1031 }
1032
1033 /* size of the Col and Row structures */
1034 #define CCOLAMD_C(n_col,ok) \
1035 ((t_mult (t_add (n_col, 1, ok), sizeof (CColamd_Col), ok) / sizeof (Int)))
1036
1037 #define CCOLAMD_R(n_row,ok) \
1038 ((t_mult (t_add (n_row, 1, ok), sizeof (CColamd_Row), ok) / sizeof (Int)))
1039
1040 /*
1041 #define CCOLAMD_RECOMMENDED(nnz, n_row, n_col) \
1042 MAX (2 * nnz, 4 * n_col) + \
1043 CCOLAMD_C (n_col) + CCOLAMD_R (n_row) + n_col + (nnz / 5) \
1044 + ((3 * n_col) + 1) + 5 * (n_col + 1) + n_row
1045 */
1046
ccolamd_need(Int nnz,Int n_row,Int n_col,int * ok)1047 static size_t ccolamd_need (Int nnz, Int n_row, Int n_col, int *ok)
1048 {
1049
1050 /* ccolamd_need, compute the following, and check for integer overflow:
1051 need = MAX (2*nnz, 4*n_col) + n_col +
1052 Col_size + Row_size +
1053 (3*n_col+1) + (5*(n_col+1)) + n_row ;
1054 */
1055 size_t s, c, r, t ;
1056
1057 /* MAX (2*nnz, 4*n_col) */
1058 s = t_mult (nnz, 2, ok) ; /* 2*nnz */
1059 t = t_mult (n_col, 4, ok) ; /* 4*n_col */
1060 s = MAX (s,t) ;
1061
1062 s = t_add (s, n_col, ok) ; /* bare minimum elbow room */
1063
1064 /* Col and Row arrays */
1065 c = CCOLAMD_C (n_col, ok) ; /* size of column structures */
1066 r = CCOLAMD_R (n_row, ok) ; /* size of row structures */
1067 s = t_add (s, c, ok) ;
1068 s = t_add (s, r, ok) ;
1069
1070 c = t_mult (n_col, 3, ok) ; /* 3*n_col + 1 */
1071 c = t_add (c, 1, ok) ;
1072 s = t_add (s, c, ok) ;
1073
1074 c = t_add (n_col, 1, ok) ; /* 5 * (n_col + 1) */
1075 c = t_mult (c, 5, ok) ;
1076 s = t_add (s, c, ok) ;
1077
1078 s = t_add (s, n_row, ok) ; /* n_row */
1079
1080 return (ok ? s : 0) ;
1081 }
1082
CCOLAMD_recommended(Int nnz,Int n_row,Int n_col)1083 PUBLIC size_t CCOLAMD_recommended /* returns recommended value of Alen. */
1084 (
1085 /* === Parameters ======================================================= */
1086
1087 Int nnz, /* number of nonzeros in A */
1088 Int n_row, /* number of rows in A */
1089 Int n_col /* number of columns in A */
1090 )
1091 {
1092 size_t s ;
1093 int ok = TRUE ;
1094 if (nnz < 0 || n_row < 0 || n_col < 0)
1095 {
1096 return (0) ;
1097 }
1098 s = ccolamd_need (nnz, n_row, n_col, &ok) ; /* bare minimum needed */
1099 s = t_add (s, nnz/5, &ok) ; /* extra elbow room */
1100 ok = ok && (s < Int_MAX) ;
1101 return (ok ? s : 0) ;
1102 }
1103
1104
1105 /* ========================================================================== */
1106 /* === ccolamd_set_defaults ================================================= */
1107 /* ========================================================================== */
1108
1109 /*
1110 * The ccolamd_set_defaults routine sets the default values of the user-
1111 * controllable parameters for ccolamd.
1112 */
1113
CCOLAMD_set_defaults(double knobs[CCOLAMD_KNOBS])1114 PUBLIC void CCOLAMD_set_defaults
1115 (
1116 /* === Parameters ======================================================= */
1117
1118 double knobs [CCOLAMD_KNOBS] /* knob array */
1119 )
1120 {
1121 /* === Local variables ================================================== */
1122
1123 Int i ;
1124
1125 if (!knobs)
1126 {
1127 return ; /* no knobs to initialize */
1128 }
1129 for (i = 0 ; i < CCOLAMD_KNOBS ; i++)
1130 {
1131 knobs [i] = 0 ;
1132 }
1133 knobs [CCOLAMD_DENSE_ROW] = 10 ;
1134 knobs [CCOLAMD_DENSE_COL] = 10 ;
1135 knobs [CCOLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
1136 knobs [CCOLAMD_LU] = FALSE ; /* default: order for Cholesky */
1137 }
1138
1139
1140 /* ========================================================================== */
1141 /* === symamd =============================================================== */
1142 /* ========================================================================== */
1143
CSYMAMD_MAIN(Int n,Int A[],Int p[],Int perm[],double knobs[CCOLAMD_KNOBS],Int stats[CCOLAMD_STATS],void * (* allocate)(size_t,size_t),void (* release)(void *),Int cmember[],Int stype)1144 PUBLIC Int CSYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */
1145 (
1146 /* === Parameters ======================================================= */
1147
1148 Int n, /* number of rows and columns of A */
1149 Int A [ ], /* row indices of A */
1150 Int p [ ], /* column pointers of A */
1151 Int perm [ ], /* output permutation, size n+1 */
1152 double knobs [CCOLAMD_KNOBS], /* parameters (uses defaults if NULL) */
1153 Int stats [CCOLAMD_STATS], /* output statistics and error codes */
1154 void * (*allocate) (size_t, size_t),/* pointer to calloc (ANSI C) or */
1155 /* mxCalloc (for MATLAB mexFunction) */
1156 void (*release) (void *), /* pointer to free (ANSI C) or */
1157 /* mxFree (for MATLAB mexFunction) */
1158 Int cmember [ ], /* constraint set */
1159 Int stype /* stype of A */
1160 )
1161 {
1162 /* === Local variables ================================================== */
1163
1164 double cknobs [CCOLAMD_KNOBS] ;
1165 double default_knobs [CCOLAMD_KNOBS] ;
1166
1167 Int *count ; /* length of each column of M, and col pointer*/
1168 Int *mark ; /* mark array for finding duplicate entries */
1169 Int *M ; /* row indices of matrix M */
1170 size_t Mlen ; /* length of M */
1171 Int n_row ; /* number of rows in M */
1172 Int nnz ; /* number of entries in A */
1173 Int i ; /* row index of A */
1174 Int j ; /* column index of A */
1175 Int k ; /* row index of M */
1176 Int mnz ; /* number of nonzeros in M */
1177 Int pp ; /* index into a column of A */
1178 Int last_row ; /* last row seen in the current column */
1179 Int length ; /* number of nonzeros in a column */
1180 Int both ; /* TRUE if ordering A+A' */
1181 Int upper ; /* TRUE if ordering triu(A)+triu(A)' */
1182 Int lower ; /* TRUE if ordering tril(A)+tril(A)' */
1183
1184 #ifndef NDEBUG
1185 ccolamd_get_debug ("csymamd") ;
1186 #endif
1187
1188 both = (stype == 0) ;
1189 upper = (stype > 0) ;
1190 lower = (stype < 0) ;
1191
1192 /* === Check the input arguments ======================================== */
1193
1194 if (!stats)
1195 {
1196 DEBUG1 (("csymamd: stats not present\n")) ;
1197 return (FALSE) ;
1198 }
1199 for (i = 0 ; i < CCOLAMD_STATS ; i++)
1200 {
1201 stats [i] = 0 ;
1202 }
1203 stats [CCOLAMD_STATUS] = CCOLAMD_OK ;
1204 stats [CCOLAMD_INFO1] = -1 ;
1205 stats [CCOLAMD_INFO2] = -1 ;
1206
1207 if (!A)
1208 {
1209 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_not_present ;
1210 DEBUG1 (("csymamd: A not present\n")) ;
1211 return (FALSE) ;
1212 }
1213
1214 if (!p) /* p is not present */
1215 {
1216 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p_not_present ;
1217 DEBUG1 (("csymamd: p not present\n")) ;
1218 return (FALSE) ;
1219 }
1220
1221 if (n < 0) /* n must be >= 0 */
1222 {
1223 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_ncol_negative ;
1224 stats [CCOLAMD_INFO1] = n ;
1225 DEBUG1 (("csymamd: n negative "ID" \n", n)) ;
1226 return (FALSE) ;
1227 }
1228
1229 nnz = p [n] ;
1230 if (nnz < 0) /* nnz must be >= 0 */
1231 {
1232 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nnz_negative ;
1233 stats [CCOLAMD_INFO1] = nnz ;
1234 DEBUG1 (("csymamd: number of entries negative "ID" \n", nnz)) ;
1235 return (FALSE) ;
1236 }
1237
1238 if (p [0] != 0)
1239 {
1240 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p0_nonzero ;
1241 stats [CCOLAMD_INFO1] = p [0] ;
1242 DEBUG1 (("csymamd: p[0] not zero "ID"\n", p [0])) ;
1243 return (FALSE) ;
1244 }
1245
1246 /* === If no knobs, set default knobs =================================== */
1247
1248 if (!knobs)
1249 {
1250 CCOLAMD_set_defaults (default_knobs) ;
1251 knobs = default_knobs ;
1252 }
1253
1254 /* === Allocate count and mark ========================================== */
1255
1256 count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1257 if (!count)
1258 {
1259 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
1260 DEBUG1 (("csymamd: allocate count (size "ID") failed\n", n+1)) ;
1261 return (FALSE) ;
1262 }
1263
1264 mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
1265 if (!mark)
1266 {
1267 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
1268 (*release) ((void *) count) ;
1269 DEBUG1 (("csymamd: allocate mark (size "ID") failed\n", n+1)) ;
1270 return (FALSE) ;
1271 }
1272
1273 /* === Compute column counts of M, check if A is valid ================== */
1274
1275 stats [CCOLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
1276
1277 for (i = 0 ; i < n ; i++)
1278 {
1279 mark [i] = -1 ;
1280 }
1281
1282 for (j = 0 ; j < n ; j++)
1283 {
1284 last_row = -1 ;
1285
1286 length = p [j+1] - p [j] ;
1287 if (length < 0)
1288 {
1289 /* column pointers must be non-decreasing */
1290 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_col_length_negative ;
1291 stats [CCOLAMD_INFO1] = j ;
1292 stats [CCOLAMD_INFO2] = length ;
1293 (*release) ((void *) count) ;
1294 (*release) ((void *) mark) ;
1295 DEBUG1 (("csymamd: col "ID" negative length "ID"\n", j, length)) ;
1296 return (FALSE) ;
1297 }
1298
1299 for (pp = p [j] ; pp < p [j+1] ; pp++)
1300 {
1301 i = A [pp] ;
1302 if (i < 0 || i >= n)
1303 {
1304 /* row index i, in column j, is out of bounds */
1305 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_row_index_out_of_bounds ;
1306 stats [CCOLAMD_INFO1] = j ;
1307 stats [CCOLAMD_INFO2] = i ;
1308 stats [CCOLAMD_INFO3] = n ;
1309 (*release) ((void *) count) ;
1310 (*release) ((void *) mark) ;
1311 DEBUG1 (("csymamd: row "ID" col "ID" out of bounds\n", i, j)) ;
1312 return (FALSE) ;
1313 }
1314
1315 if (i <= last_row || mark [i] == j)
1316 {
1317 /* row index is unsorted or repeated (or both), thus col */
1318 /* is jumbled. This is a notice, not an error condition. */
1319 stats [CCOLAMD_STATUS] = CCOLAMD_OK_BUT_JUMBLED ;
1320 stats [CCOLAMD_INFO1] = j ;
1321 stats [CCOLAMD_INFO2] = i ;
1322 (stats [CCOLAMD_INFO3]) ++ ;
1323 DEBUG1 (("csymamd: row "ID" col "ID" unsorted/dupl.\n", i, j)) ;
1324 }
1325
1326 if (mark [i] != j)
1327 {
1328 if ((both && i != j) || (lower && i > j) || (upper && i < j))
1329 {
1330 /* row k of M will contain column indices i and j */
1331 count [i]++ ;
1332 count [j]++ ;
1333 }
1334 }
1335
1336 /* mark the row as having been seen in this column */
1337 mark [i] = j ;
1338
1339 last_row = i ;
1340 }
1341 }
1342
1343 /* === Compute column pointers of M ===================================== */
1344
1345 /* use output permutation, perm, for column pointers of M */
1346 perm [0] = 0 ;
1347 for (j = 1 ; j <= n ; j++)
1348 {
1349 perm [j] = perm [j-1] + count [j-1] ;
1350 }
1351 for (j = 0 ; j < n ; j++)
1352 {
1353 count [j] = perm [j] ;
1354 }
1355
1356 /* === Construct M ====================================================== */
1357
1358 mnz = perm [n] ;
1359 n_row = mnz / 2 ;
1360 Mlen = CCOLAMD_recommended (mnz, n_row, n) ;
1361 M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
1362 DEBUG1 (("csymamd: M is "ID"-by-"ID" with "ID" entries, Mlen = %g\n",
1363 n_row, n, mnz, (double) Mlen)) ;
1364
1365 if (!M)
1366 {
1367 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
1368 (*release) ((void *) count) ;
1369 (*release) ((void *) mark) ;
1370 DEBUG1 (("csymamd: allocate M (size %g) failed\n", (double) Mlen)) ;
1371 return (FALSE) ;
1372 }
1373
1374 k = 0 ;
1375
1376 if (stats [CCOLAMD_STATUS] == CCOLAMD_OK)
1377 {
1378 /* Matrix is OK */
1379 for (j = 0 ; j < n ; j++)
1380 {
1381 ASSERT (p [j+1] - p [j] >= 0) ;
1382 for (pp = p [j] ; pp < p [j+1] ; pp++)
1383 {
1384 i = A [pp] ;
1385 ASSERT (i >= 0 && i < n) ;
1386 if ((both && i != j) || (lower && i > j) || (upper && i < j))
1387 {
1388 /* row k of M contains column indices i and j */
1389 M [count [i]++] = k ;
1390 M [count [j]++] = k ;
1391 k++ ;
1392 }
1393 }
1394 }
1395 }
1396 else
1397 {
1398 /* Matrix is jumbled. Do not add duplicates to M. Unsorted cols OK. */
1399 DEBUG1 (("csymamd: Duplicates in A.\n")) ;
1400 for (i = 0 ; i < n ; i++)
1401 {
1402 mark [i] = -1 ;
1403 }
1404 for (j = 0 ; j < n ; j++)
1405 {
1406 ASSERT (p [j+1] - p [j] >= 0) ;
1407 for (pp = p [j] ; pp < p [j+1] ; pp++)
1408 {
1409 i = A [pp] ;
1410 ASSERT (i >= 0 && i < n) ;
1411 if (mark [i] != j)
1412 {
1413 if ((both && i != j) || (lower && i > j) || (upper && i<j))
1414 {
1415 /* row k of M contains column indices i and j */
1416 M [count [i]++] = k ;
1417 M [count [j]++] = k ;
1418 k++ ;
1419 mark [i] = j ;
1420 }
1421 }
1422 }
1423 }
1424 }
1425
1426 /* count and mark no longer needed */
1427 (*release) ((void *) mark) ;
1428 (*release) ((void *) count) ;
1429 ASSERT (k == n_row) ;
1430
1431 /* === Adjust the knobs for M =========================================== */
1432
1433 for (i = 0 ; i < CCOLAMD_KNOBS ; i++)
1434 {
1435 cknobs [i] = knobs [i] ;
1436 }
1437
1438 /* there are no dense rows in M */
1439 cknobs [CCOLAMD_DENSE_ROW] = -1 ;
1440 cknobs [CCOLAMD_DENSE_COL] = knobs [CCOLAMD_DENSE_ROW] ;
1441
1442 /* ensure CCSYMAMD orders for Cholesky, not LU */
1443 cknobs [CCOLAMD_LU] = FALSE ;
1444
1445 /* === Order the columns of M =========================================== */
1446
1447 (void) CCOLAMD_2 (n_row, n, (Int) Mlen, M, perm, cknobs, stats,
1448 (Int *) NULL, (Int *) NULL, (Int *) NULL, (Int *) NULL,
1449 (Int *) NULL, (Int *) NULL, (Int *) NULL, cmember) ;
1450
1451 /* === adjust statistics ================================================ */
1452
1453 /* a dense column in ccolamd means a dense row and col in csymamd */
1454 stats [CCOLAMD_DENSE_ROW] = stats [CCOLAMD_DENSE_COL] ;
1455
1456 /* === Free M =========================================================== */
1457
1458 (*release) ((void *) M) ;
1459 DEBUG1 (("csymamd: done.\n")) ;
1460 return (TRUE) ;
1461 }
1462
1463
1464 /* ========================================================================== */
1465 /* === ccolamd ============================================================== */
1466 /* ========================================================================== */
1467
1468 /*
1469 * The colamd routine computes a column ordering Q of a sparse matrix
1470 * A such that the LU factorization P(AQ) = LU remains sparse, where P is
1471 * selected via partial pivoting. The routine can also be viewed as
1472 * providing a permutation Q such that the Cholesky factorization
1473 * (AQ)'(AQ) = LL' remains sparse.
1474 */
1475
CCOLAMD_MAIN(Int n_row,Int n_col,Int Alen,Int A[],Int p[],double knobs[CCOLAMD_KNOBS],Int stats[CCOLAMD_STATS],Int cmember[])1476 PUBLIC Int CCOLAMD_MAIN
1477 (
1478 /* === Parameters ======================================================= */
1479
1480 Int n_row, /* number of rows in A */
1481 Int n_col, /* number of columns in A */
1482 Int Alen, /* length of A */
1483 Int A [ ], /* row indices of A */
1484 Int p [ ], /* pointers to columns in A */
1485 double knobs [CCOLAMD_KNOBS],/* parameters (uses defaults if NULL) */
1486 Int stats [CCOLAMD_STATS], /* output statistics and error codes */
1487 Int cmember [ ] /* constraint set of A */
1488 )
1489 {
1490 return (CCOLAMD_2 (n_row, n_col, Alen, A, p, knobs, stats,
1491 (Int *) NULL, (Int *) NULL, (Int *) NULL, (Int *) NULL,
1492 (Int *) NULL, (Int *) NULL, (Int *) NULL, cmember)) ;
1493 }
1494
1495
1496 /* ========================================================================== */
1497 /* === ccolamd2 ============================================================= */
1498 /* ========================================================================== */
1499
1500 /* Identical to ccolamd, except that additional information about each frontal
1501 * matrix is returned to the caller. Not intended to be directly called by
1502 * the user.
1503 */
1504
CCOLAMD_2(Int n_row,Int n_col,Int Alen,Int A[],Int p[],double knobs[CCOLAMD_KNOBS],Int stats[CCOLAMD_STATS],Int Front_npivcol[],Int Front_nrows[],Int Front_ncols[],Int Front_parent[],Int Front_cols[],Int * p_nfr,Int InFront[],Int cmember[])1505 PUBLIC Int CCOLAMD_2 /* returns TRUE if successful, FALSE otherwise */
1506 (
1507 /* === Parameters ======================================================= */
1508
1509 Int n_row, /* number of rows in A */
1510 Int n_col, /* number of columns in A */
1511 Int Alen, /* length of A */
1512 Int A [ ], /* row indices of A */
1513 Int p [ ], /* pointers to columns in A */
1514 double knobs [CCOLAMD_KNOBS],/* parameters (uses defaults if NULL) */
1515 Int stats [CCOLAMD_STATS], /* output statistics and error codes */
1516
1517 /* each Front array is of size n_col+1. */
1518 Int Front_npivcol [ ], /* # pivot cols in each front */
1519 Int Front_nrows [ ], /* # of rows in each front (incl. pivot rows) */
1520 Int Front_ncols [ ], /* # of cols in each front (incl. pivot cols) */
1521 Int Front_parent [ ], /* parent of each front */
1522 Int Front_cols [ ], /* link list of pivot columns for each front */
1523 Int *p_nfr, /* total number of frontal matrices */
1524 Int InFront [ ], /* InFront [row] = f if the original row was
1525 * absorbed into front f. EMPTY if the row was
1526 * empty, dense, or not absorbed. This array
1527 * has size n_row+1 */
1528 Int cmember [ ] /* constraint set of A */
1529 )
1530 {
1531 /* === Local variables ================================================== */
1532
1533 Int i ; /* loop index */
1534 Int nnz ; /* nonzeros in A */
1535 size_t Row_size ; /* size of Row [ ], in integers */
1536 size_t Col_size ; /* size of Col [ ], in integers */
1537 size_t need ; /* minimum required length of A */
1538 CColamd_Row *Row ; /* pointer into A of Row [0..n_row] array */
1539 CColamd_Col *Col ; /* pointer into A of Col [0..n_col] array */
1540 Int n_col2 ; /* number of non-dense, non-empty columns */
1541 Int n_row2 ; /* number of non-dense, non-empty rows */
1542 Int ngarbage ; /* number of garbage collections performed */
1543 Int max_deg ; /* maximum row degree */
1544 double default_knobs [CCOLAMD_KNOBS] ; /* default knobs array */
1545
1546 Int n_cset ; /* number of constraint sets */
1547 Int *cset ; /* cset of A */
1548 Int *cset_start ; /* pointer into cset */
1549 Int *temp_cstart ; /* temp pointer to start of cset */
1550 Int *csize ; /* temp pointer to cset size */
1551 Int ap ; /* column index */
1552 Int order_for_lu ; /* TRUE: order for LU, FALSE: for Cholesky */
1553
1554 Int ndense_row, nempty_row, parent, ndense_col,
1555 nempty_col, k, col, nfr, *Front_child, *Front_sibling, *Front_stack,
1556 *Front_order, *Front_size ;
1557 Int nnewlyempty_col, nnewlyempty_row ;
1558 Int aggressive ;
1559 Int row ;
1560 Int *dead_cols ;
1561 Int set1 ;
1562 Int set2 ;
1563 Int cs ;
1564
1565 int ok ;
1566
1567 #ifndef NDEBUG
1568 ccolamd_get_debug ("ccolamd") ;
1569 #endif
1570
1571 /* === Check the input arguments ======================================== */
1572
1573 if (!stats)
1574 {
1575 DEBUG1 (("ccolamd: stats not present\n")) ;
1576 return (FALSE) ;
1577 }
1578 for (i = 0 ; i < CCOLAMD_STATS ; i++)
1579 {
1580 stats [i] = 0 ;
1581 }
1582 stats [CCOLAMD_STATUS] = CCOLAMD_OK ;
1583 stats [CCOLAMD_INFO1] = -1 ;
1584 stats [CCOLAMD_INFO2] = -1 ;
1585
1586 if (!A) /* A is not present */
1587 {
1588 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_not_present ;
1589 DEBUG1 (("ccolamd: A not present\n")) ;
1590 return (FALSE) ;
1591 }
1592
1593 if (!p) /* p is not present */
1594 {
1595 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p_not_present ;
1596 DEBUG1 (("ccolamd: p not present\n")) ;
1597 return (FALSE) ;
1598 }
1599
1600 if (n_row < 0) /* n_row must be >= 0 */
1601 {
1602 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nrow_negative ;
1603 stats [CCOLAMD_INFO1] = n_row ;
1604 DEBUG1 (("ccolamd: nrow negative "ID"\n", n_row)) ;
1605 return (FALSE) ;
1606 }
1607
1608 if (n_col < 0) /* n_col must be >= 0 */
1609 {
1610 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_ncol_negative ;
1611 stats [CCOLAMD_INFO1] = n_col ;
1612 DEBUG1 (("ccolamd: ncol negative "ID"\n", n_col)) ;
1613 return (FALSE) ;
1614 }
1615
1616 nnz = p [n_col] ;
1617 if (nnz < 0) /* nnz must be >= 0 */
1618 {
1619 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nnz_negative ;
1620 stats [CCOLAMD_INFO1] = nnz ;
1621 DEBUG1 (("ccolamd: number of entries negative "ID"\n", nnz)) ;
1622 return (FALSE) ;
1623 }
1624
1625 if (p [0] != 0)
1626 {
1627 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p0_nonzero ;
1628 stats [CCOLAMD_INFO1] = p [0] ;
1629 DEBUG1 (("ccolamd: p[0] not zero "ID"\n", p [0])) ;
1630 return (FALSE) ;
1631 }
1632
1633 /* === If no knobs, set default knobs =================================== */
1634
1635 if (!knobs)
1636 {
1637 CCOLAMD_set_defaults (default_knobs) ;
1638 knobs = default_knobs ;
1639 }
1640
1641 aggressive = (knobs [CCOLAMD_AGGRESSIVE] != FALSE) ;
1642 order_for_lu = (knobs [CCOLAMD_LU] != FALSE) ;
1643
1644 /* === Allocate workspace from array A ================================== */
1645
1646 ok = TRUE ;
1647 Col_size = CCOLAMD_C (n_col, &ok) ;
1648 Row_size = CCOLAMD_R (n_row, &ok) ;
1649
1650 /* min size of A is 2nnz+ncol. cset and cset_start are of size 2ncol+1 */
1651 /* Each of the 5 fronts is of size n_col + 1. InFront is of size nrow. */
1652
1653 /*
1654 need = MAX (2*nnz, 4*n_col) + n_col +
1655 Col_size + Row_size +
1656 (3*n_col+1) + (5*(n_col+1)) + n_row ;
1657 */
1658 need = ccolamd_need (nnz, n_row, n_col, &ok) ;
1659
1660 if (!ok || need > (size_t) Alen || need > Int_MAX)
1661 {
1662 /* not enough space in array A to perform the ordering */
1663 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_too_small ;
1664 stats [CCOLAMD_INFO1] = need ;
1665 stats [CCOLAMD_INFO2] = Alen ;
1666 DEBUG1 (("ccolamd: Need Alen >= "ID", given "ID"\n", need, Alen)) ;
1667 return (FALSE) ;
1668 }
1669
1670 /* since integer overflow has been check, the following cannot overflow: */
1671 Alen -= Col_size + Row_size + (3*n_col + 1) + 5*(n_col+1) + n_row ;
1672
1673 /* Size of A is now Alen >= MAX (2*nnz, 4*n_col) + n_col. The ordering
1674 * requires Alen >= 2*nnz + n_col, and the postorder requires
1675 * Alen >= 5*n_col. */
1676
1677 ap = Alen ;
1678
1679 /* Front array workspace: 5*(n_col+1) + n_row */
1680 if (!Front_npivcol || !Front_nrows || !Front_ncols || !Front_parent ||
1681 !Front_cols || !Front_cols || !InFront)
1682 {
1683 Front_npivcol = &A [ap] ; ap += (n_col + 1) ;
1684 Front_nrows = &A [ap] ; ap += (n_col + 1) ;
1685 Front_ncols = &A [ap] ; ap += (n_col + 1) ;
1686 Front_parent = &A [ap] ; ap += (n_col + 1) ;
1687 Front_cols = &A [ap] ; ap += (n_col + 1) ;
1688 InFront = &A [ap] ; ap += (n_row) ;
1689 }
1690 else
1691 {
1692 /* Fronts are present. Leave the additional space as elbow room. */
1693 ap += 5*(n_col+1) + n_row ;
1694 ap = Alen ;
1695 }
1696
1697 /* Workspace for cset management: 3*n_col+1 */
1698 /* cset_start is of size n_col + 1 */
1699 cset_start = &A [ap] ;
1700 ap += n_col + 1 ;
1701
1702 /* dead_col is of size n_col */
1703 dead_cols = &A [ap] ;
1704 ap += n_col ;
1705
1706 /* cset is of size n_col */
1707 cset = &A [ap] ;
1708 ap += n_col ;
1709
1710 /* Col is of size Col_size. The space is shared by temp_cstart and csize */
1711 Col = (CColamd_Col *) &A [ap] ;
1712 temp_cstart = (Int *) Col ; /* [ temp_cstart is of size n_col+1 */
1713 csize = temp_cstart + (n_col+1) ; /* csize is of size n_col+1 */
1714 ap += Col_size ;
1715 ASSERT (Col_size >= 2*n_col+1) ;
1716
1717 /* Row is of size Row_size */
1718 Row = (CColamd_Row *) &A [ap] ;
1719 ap += Row_size ;
1720
1721 /* Initialize csize & dead_cols to zero */
1722 for (i = 0 ; i < n_col ; i++)
1723 {
1724 csize [i] = 0 ;
1725 dead_cols [i] = 0 ;
1726 }
1727
1728 /* === Construct the constraint set ===================================== */
1729
1730 if (n_col == 0)
1731 {
1732 n_cset = 0 ;
1733 }
1734 else if (cmember == (Int *) NULL)
1735 {
1736 /* no constraint set; all columns belong to set zero */
1737 n_cset = 1 ;
1738 csize [0] = n_col ;
1739 DEBUG1 (("no cmember present\n")) ;
1740 }
1741 else
1742 {
1743 n_cset = 0 ;
1744 for (i = 0 ; i < n_col ; i++)
1745 {
1746 if (cmember [i] < 0 || cmember [i] > n_col)
1747 {
1748 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_invalid_cmember ;
1749 DEBUG1 (("ccolamd: malformed cmember \n")) ;
1750 return (FALSE) ;
1751 }
1752 n_cset = MAX (n_cset, cmember [i]) ;
1753 csize [cmember [i]]++ ;
1754 }
1755 /* cset is zero based */
1756 n_cset++ ;
1757 }
1758
1759 ASSERT ((n_cset >= 0) && (n_cset <= n_col)) ;
1760
1761 cset_start [0] = temp_cstart [0] = 0 ;
1762 for (i = 1 ; i <= n_cset ; i++)
1763 {
1764 cset_start [i] = cset_start [i-1] + csize [i-1] ;
1765 DEBUG4 ((" cset_start ["ID"] = "ID" \n", i , cset_start [i])) ;
1766 temp_cstart [i] = cset_start [i] ;
1767 }
1768
1769 /* do in reverse order to encourage natural tie-breaking */
1770 if (cmember == (Int *) NULL)
1771 {
1772 for (i = n_col-1 ; i >= 0 ; i--)
1773 {
1774 cset [temp_cstart [0]++] = i ;
1775 }
1776 }
1777 else
1778 {
1779 for (i = n_col-1 ; i >= 0 ; i--)
1780 {
1781 cset [temp_cstart [cmember [i]]++] = i ;
1782 }
1783 }
1784
1785 /* ] temp_cstart and csize are no longer used */
1786
1787 /* === Construct the row and column data structures ===================== */
1788
1789 if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
1790 {
1791 /* input matrix is invalid */
1792 DEBUG1 (("ccolamd: Matrix invalid\n")) ;
1793 return (FALSE) ;
1794 }
1795
1796 /* === Initialize front info ============================================ */
1797
1798 for (col = 0 ; col < n_col ; col++)
1799 {
1800 Front_npivcol [col] = 0 ;
1801 Front_nrows [col] = 0 ;
1802 Front_ncols [col] = 0 ;
1803 Front_parent [col] = EMPTY ;
1804 Front_cols [col] = EMPTY ;
1805 }
1806
1807 /* === Initialize scores, kill dense rows/columns ======================= */
1808
1809 init_scoring (n_row, n_col, Row, Col, A, p, knobs,
1810 &n_row2, &n_col2, &max_deg, cmember, n_cset, cset_start, dead_cols,
1811 &ndense_row, &nempty_row, &nnewlyempty_row,
1812 &ndense_col, &nempty_col, &nnewlyempty_col) ;
1813
1814 ASSERT (n_row2 == n_row - nempty_row - nnewlyempty_row - ndense_row) ;
1815 ASSERT (n_col2 == n_col - nempty_col - nnewlyempty_col - ndense_col) ;
1816 DEBUG1 (("# dense rows "ID" cols "ID"\n", ndense_row, ndense_col)) ;
1817
1818 /* === Order the supercolumns =========================================== */
1819
1820 ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
1821 #ifndef NDEBUG
1822 n_col2,
1823 #endif
1824 max_deg, 2*nnz, cset, cset_start,
1825 #ifndef NDEBUG
1826 n_cset,
1827 #endif
1828 cmember, Front_npivcol, Front_nrows, Front_ncols, Front_parent,
1829 Front_cols, &nfr, aggressive, InFront, order_for_lu) ;
1830
1831 ASSERT (Alen >= 5*n_col) ;
1832
1833 /* === Postorder ======================================================== */
1834
1835 /* A is no longer needed, so use A [0..5*nfr-1] as workspace [ [ */
1836 /* This step requires Alen >= 5*n_col */
1837 Front_child = A ;
1838 Front_sibling = Front_child + nfr ;
1839 Front_stack = Front_sibling + nfr ;
1840 Front_order = Front_stack + nfr ;
1841 Front_size = Front_order + nfr ;
1842
1843 CCOLAMD_fsize (nfr, Front_size, Front_nrows, Front_ncols,
1844 Front_parent, Front_npivcol) ;
1845
1846 CCOLAMD_postorder (nfr, Front_parent, Front_npivcol, Front_size,
1847 Front_order, Front_child, Front_sibling, Front_stack, Front_cols,
1848 cmember) ;
1849
1850 /* Front_size, Front_stack, Front_child, Front_sibling no longer needed ] */
1851
1852 /* use A [0..nfr-1] as workspace */
1853 CCOLAMD_apply_order (Front_npivcol, Front_order, A, nfr, nfr) ;
1854 CCOLAMD_apply_order (Front_nrows, Front_order, A, nfr, nfr) ;
1855 CCOLAMD_apply_order (Front_ncols, Front_order, A, nfr, nfr) ;
1856 CCOLAMD_apply_order (Front_parent, Front_order, A, nfr, nfr) ;
1857 CCOLAMD_apply_order (Front_cols, Front_order, A, nfr, nfr) ;
1858
1859 /* fix the parent to refer to the new numbering */
1860 for (i = 0 ; i < nfr ; i++)
1861 {
1862 parent = Front_parent [i] ;
1863 if (parent != EMPTY)
1864 {
1865 Front_parent [i] = Front_order [parent] ;
1866 }
1867 }
1868
1869 /* fix InFront to refer to the new numbering */
1870 for (row = 0 ; row < n_row ; row++)
1871 {
1872 i = InFront [row] ;
1873 ASSERT (i >= EMPTY && i < nfr) ;
1874 if (i != EMPTY)
1875 {
1876 InFront [row] = Front_order [i] ;
1877 }
1878 }
1879
1880 /* Front_order longer needed ] */
1881
1882 /* === Order the columns in the fronts ================================== */
1883
1884 /* use A [0..n_col-1] as inverse permutation */
1885 for (i = 0 ; i < n_col ; i++)
1886 {
1887 A [i] = EMPTY ;
1888 }
1889
1890 k = 0 ;
1891 set1 = 0 ;
1892 for (i = 0 ; i < nfr ; i++)
1893 {
1894 ASSERT (Front_npivcol [i] > 0) ;
1895
1896 set2 = CMEMBER (Front_cols [i]) ;
1897 while (set1 < set2)
1898 {
1899 k += dead_cols [set1] ;
1900 DEBUG3 (("Skip null/dense columns of set "ID"\n",set1)) ;
1901 set1++ ;
1902 }
1903 set1 = set2 ;
1904
1905 for (col = Front_cols [i] ; col != EMPTY ; col = Col [col].nextcol)
1906 {
1907 ASSERT (col >= 0 && col < n_col) ;
1908 DEBUG1 (("ccolamd output ordering: k "ID" col "ID"\n", k, col)) ;
1909 p [k] = col ;
1910 ASSERT (A [col] == EMPTY) ;
1911
1912 cs = CMEMBER (col) ;
1913 ASSERT (k >= cset_start [cs] && k < cset_start [cs+1]) ;
1914
1915 A [col] = k ;
1916 k++ ;
1917 }
1918 }
1919
1920 /* === Order the "dense" and null columns =============================== */
1921
1922 if (n_col2 < n_col)
1923 {
1924 for (col = 0 ; col < n_col ; col++)
1925 {
1926 if (A [col] == EMPTY)
1927 {
1928 k = Col [col].shared2.order ;
1929 cs = CMEMBER (col) ;
1930 #ifndef NDEBUG
1931 dead_cols [cs]-- ;
1932 #endif
1933 ASSERT (k >= cset_start [cs] && k < cset_start [cs+1]) ;
1934 DEBUG1 (("ccolamd output ordering: k "ID" col "ID
1935 " (dense or null col)\n", k, col)) ;
1936 p [k] = col ;
1937 A [col] = k ;
1938 }
1939 }
1940 }
1941
1942 #ifndef NDEBUG
1943 for (i = 0 ; i < n_cset ; i++)
1944 {
1945 ASSERT (dead_cols [i] == 0) ;
1946 }
1947 #endif
1948
1949 /* === Return statistics in stats ======================================= */
1950
1951 stats [CCOLAMD_DENSE_ROW] = ndense_row ;
1952 stats [CCOLAMD_EMPTY_ROW] = nempty_row ; /* fixed in 2.7.3 */
1953 stats [CCOLAMD_NEWLY_EMPTY_ROW] = nnewlyempty_row ;
1954 stats [CCOLAMD_DENSE_COL] = ndense_col ;
1955 stats [CCOLAMD_EMPTY_COL] = nempty_col ;
1956 stats [CCOLAMD_NEWLY_EMPTY_COL] = nnewlyempty_col ;
1957 ASSERT (ndense_col + nempty_col + nnewlyempty_col == n_col - n_col2) ;
1958 if (p_nfr)
1959 {
1960 *p_nfr = nfr ;
1961 }
1962 stats [CCOLAMD_DEFRAG_COUNT] = ngarbage ;
1963 DEBUG1 (("ccolamd: done.\n")) ;
1964 return (TRUE) ;
1965 }
1966
1967
1968 /* ========================================================================== */
1969 /* === colamd_report ======================================================== */
1970 /* ========================================================================== */
1971
CCOLAMD_report(Int stats[CCOLAMD_STATS])1972 PUBLIC void CCOLAMD_report
1973 (
1974 Int stats [CCOLAMD_STATS]
1975 )
1976 {
1977 print_report ("ccolamd", stats) ;
1978 }
1979
1980
1981 /* ========================================================================== */
1982 /* === symamd_report ======================================================== */
1983 /* ========================================================================== */
1984
CSYMAMD_report(Int stats[CCOLAMD_STATS])1985 PUBLIC void CSYMAMD_report
1986 (
1987 Int stats [CCOLAMD_STATS]
1988 )
1989 {
1990 print_report ("csymamd", stats) ;
1991 }
1992
1993
1994 /* ========================================================================== */
1995 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
1996 /* ========================================================================== */
1997
1998 /* There are no user-callable routines beyond this point in the file */
1999
2000
2001 /* ========================================================================== */
2002 /* === init_rows_cols ======================================================= */
2003 /* ========================================================================== */
2004
2005 /*
2006 Takes the column form of the matrix in A and creates the row form of the
2007 matrix. Also, row and column attributes are stored in the Col and Row
2008 structs. If the columns are un-sorted or contain duplicate row indices,
2009 this routine will also sort and remove duplicate row indices from the
2010 column form of the matrix. Returns FALSE if the matrix is invalid,
2011 TRUE otherwise. Not user-callable.
2012 */
2013
init_rows_cols(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int p[],Int stats[CCOLAMD_STATS])2014 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
2015 (
2016 /* === Parameters ======================================================= */
2017
2018 Int n_row, /* number of rows of A */
2019 Int n_col, /* number of columns of A */
2020 CColamd_Row Row [ ], /* of size n_row+1 */
2021 CColamd_Col Col [ ], /* of size n_col+1 */
2022 Int A [ ], /* row indices of A, of size Alen */
2023 Int p [ ], /* pointers to columns in A, of size n_col+1 */
2024 Int stats [CCOLAMD_STATS] /* colamd statistics */
2025 )
2026 {
2027 /* === Local variables ================================================== */
2028
2029 Int col ; /* a column index */
2030 Int row ; /* a row index */
2031 Int *cp ; /* a column pointer */
2032 Int *cp_end ; /* a pointer to the end of a column */
2033 Int *rp ; /* a row pointer */
2034 Int *rp_end ; /* a pointer to the end of a row */
2035 Int last_row ; /* previous row */
2036
2037 /* === Initialize columns, and check column pointers ==================== */
2038
2039 for (col = 0 ; col < n_col ; col++)
2040 {
2041 Col [col].start = p [col] ;
2042 Col [col].length = p [col+1] - p [col] ;
2043
2044 if (Col [col].length < 0)
2045 {
2046 /* column pointers must be non-decreasing */
2047 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_col_length_negative ;
2048 stats [CCOLAMD_INFO1] = col ;
2049 stats [CCOLAMD_INFO2] = Col [col].length ;
2050 DEBUG1 (("ccolamd: col "ID" length "ID" < 0\n",
2051 col, Col [col].length)) ;
2052 return (FALSE) ;
2053 }
2054
2055 Col [col].shared1.thickness = 1 ;
2056 Col [col].shared2.score = 0 ;
2057 Col [col].shared3.prev = EMPTY ;
2058 Col [col].shared4.degree_next = EMPTY ;
2059 Col [col].nextcol = EMPTY ;
2060 Col [col].lastcol = col ;
2061 }
2062
2063 /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
2064
2065 /* === Scan columns, compute row degrees, and check row indices ========= */
2066
2067 stats [CCOLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
2068
2069 for (row = 0 ; row < n_row ; row++)
2070 {
2071 Row [row].length = 0 ;
2072 Row [row].shared2.mark = -1 ;
2073 Row [row].thickness = 1 ;
2074 Row [row].front = EMPTY ;
2075 }
2076
2077 for (col = 0 ; col < n_col ; col++)
2078 {
2079 DEBUG1 (("\nCcolamd input column "ID":\n", col)) ;
2080 last_row = -1 ;
2081
2082 cp = &A [p [col]] ;
2083 cp_end = &A [p [col+1]] ;
2084
2085 while (cp < cp_end)
2086 {
2087 row = *cp++ ;
2088 DEBUG1 (("row: "ID"\n", row)) ;
2089
2090 /* make sure row indices within range */
2091 if (row < 0 || row >= n_row)
2092 {
2093 stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_row_index_out_of_bounds ;
2094 stats [CCOLAMD_INFO1] = col ;
2095 stats [CCOLAMD_INFO2] = row ;
2096 stats [CCOLAMD_INFO3] = n_row ;
2097 DEBUG1 (("row "ID" col "ID" out of bounds\n", row, col)) ;
2098 return (FALSE) ;
2099 }
2100
2101 if (row <= last_row || Row [row].shared2.mark == col)
2102 {
2103 /* row index are unsorted or repeated (or both), thus col */
2104 /* is jumbled. This is a notice, not an error condition. */
2105 stats [CCOLAMD_STATUS] = CCOLAMD_OK_BUT_JUMBLED ;
2106 stats [CCOLAMD_INFO1] = col ;
2107 stats [CCOLAMD_INFO2] = row ;
2108 (stats [CCOLAMD_INFO3]) ++ ;
2109 DEBUG1 (("row "ID" col "ID" unsorted/duplicate\n", row, col)) ;
2110 }
2111
2112 if (Row [row].shared2.mark != col)
2113 {
2114 Row [row].length++ ;
2115 }
2116 else
2117 {
2118 /* this is a repeated entry in the column, */
2119 /* it will be removed */
2120 Col [col].length-- ;
2121 }
2122
2123 /* mark the row as having been seen in this column */
2124 Row [row].shared2.mark = col ;
2125
2126 last_row = row ;
2127 }
2128 }
2129
2130 /* === Compute row pointers ============================================= */
2131
2132 /* row form of the matrix starts directly after the column */
2133 /* form of matrix in A */
2134 Row [0].start = p [n_col] ;
2135 Row [0].shared1.p = Row [0].start ;
2136 Row [0].shared2.mark = -1 ;
2137 for (row = 1 ; row < n_row ; row++)
2138 {
2139 Row [row].start = Row [row-1].start + Row [row-1].length ;
2140 Row [row].shared1.p = Row [row].start ;
2141 Row [row].shared2.mark = -1 ;
2142 }
2143
2144 /* === Create row form ================================================== */
2145
2146 if (stats [CCOLAMD_STATUS] == CCOLAMD_OK_BUT_JUMBLED)
2147 {
2148 /* if cols jumbled, watch for repeated row indices */
2149 for (col = 0 ; col < n_col ; col++)
2150 {
2151 cp = &A [p [col]] ;
2152 cp_end = &A [p [col+1]] ;
2153 while (cp < cp_end)
2154 {
2155 row = *cp++ ;
2156 if (Row [row].shared2.mark != col)
2157 {
2158 A [(Row [row].shared1.p)++] = col ;
2159 Row [row].shared2.mark = col ;
2160 }
2161 }
2162 }
2163 }
2164 else
2165 {
2166 /* if cols not jumbled, we don't need the mark (this is faster) */
2167 for (col = 0 ; col < n_col ; col++)
2168 {
2169 cp = &A [p [col]] ;
2170 cp_end = &A [p [col+1]] ;
2171 while (cp < cp_end)
2172 {
2173 A [(Row [*cp++].shared1.p)++] = col ;
2174 }
2175 }
2176 }
2177
2178 /* === Clear the row marks and set row degrees ========================== */
2179
2180 for (row = 0 ; row < n_row ; row++)
2181 {
2182 Row [row].shared2.mark = 0 ;
2183 Row [row].shared1.degree = Row [row].length ;
2184 }
2185
2186 /* === See if we need to re-create columns ============================== */
2187
2188 if (stats [CCOLAMD_STATUS] == CCOLAMD_OK_BUT_JUMBLED)
2189 {
2190 DEBUG1 (("ccolamd: reconstructing column form, matrix jumbled\n")) ;
2191
2192 #ifndef NDEBUG
2193 /* make sure column lengths are correct */
2194 for (col = 0 ; col < n_col ; col++)
2195 {
2196 p [col] = Col [col].length ;
2197 }
2198 for (row = 0 ; row < n_row ; row++)
2199 {
2200 rp = &A [Row [row].start] ;
2201 rp_end = rp + Row [row].length ;
2202 while (rp < rp_end)
2203 {
2204 p [*rp++]-- ;
2205 }
2206 }
2207 for (col = 0 ; col < n_col ; col++)
2208 {
2209 ASSERT (p [col] == 0) ;
2210 }
2211 /* now p is all zero (different than when debugging is turned off) */
2212 #endif
2213
2214 /* === Compute col pointers ========================================= */
2215
2216 /* col form of the matrix starts at A [0]. */
2217 /* Note, we may have a gap between the col form and the row */
2218 /* form if there were duplicate entries, if so, it will be */
2219 /* removed upon the first garbage collection */
2220 Col [0].start = 0 ;
2221 p [0] = Col [0].start ;
2222 for (col = 1 ; col < n_col ; col++)
2223 {
2224 /* note that the lengths here are for pruned columns, i.e. */
2225 /* no duplicate row indices will exist for these columns */
2226 Col [col].start = Col [col-1].start + Col [col-1].length ;
2227 p [col] = Col [col].start ;
2228 }
2229
2230 /* === Re-create col form =========================================== */
2231
2232 for (row = 0 ; row < n_row ; row++)
2233 {
2234 rp = &A [Row [row].start] ;
2235 rp_end = rp + Row [row].length ;
2236 while (rp < rp_end)
2237 {
2238 A [(p [*rp++])++] = row ;
2239 }
2240 }
2241 }
2242
2243 /* === Done. Matrix is not (or no longer) jumbled ====================== */
2244
2245
2246 return (TRUE) ;
2247 }
2248
2249
2250 /* ========================================================================== */
2251 /* === init_scoring ========================================================= */
2252 /* ========================================================================== */
2253
2254 /*
2255 Kills dense or empty columns and rows, calculates an initial score for
2256 each column, and places all columns in the degree lists. Not user-callable.
2257 */
2258
init_scoring(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int head[],double knobs[CCOLAMD_KNOBS],Int * p_n_row2,Int * p_n_col2,Int * p_max_deg,Int cmember[],Int n_cset,Int cset_start[],Int dead_cols[],Int * p_ndense_row,Int * p_nempty_row,Int * p_nnewlyempty_row,Int * p_ndense_col,Int * p_nempty_col,Int * p_nnewlyempty_col)2259 PRIVATE void init_scoring
2260 (
2261 /* === Parameters ======================================================= */
2262
2263 Int n_row, /* number of rows of A */
2264 Int n_col, /* number of columns of A */
2265 CColamd_Row Row [ ], /* of size n_row+1 */
2266 CColamd_Col Col [ ], /* of size n_col+1 */
2267 Int A [ ], /* column form and row form of A */
2268 Int head [ ], /* of size n_col+1 */
2269 double knobs [CCOLAMD_KNOBS],/* parameters */
2270 Int *p_n_row2, /* number of non-dense, non-empty rows */
2271 Int *p_n_col2, /* number of non-dense, non-empty columns */
2272 Int *p_max_deg, /* maximum row degree */
2273 Int cmember [ ],
2274 Int n_cset,
2275 Int cset_start [ ],
2276 Int dead_cols [ ],
2277 Int *p_ndense_row, /* number of dense rows */
2278 Int *p_nempty_row, /* number of original empty rows */
2279 Int *p_nnewlyempty_row, /* number of newly empty rows */
2280 Int *p_ndense_col, /* number of dense cols (excl "empty" cols) */
2281 Int *p_nempty_col, /* number of original empty cols */
2282 Int *p_nnewlyempty_col /* number of newly empty cols */
2283 )
2284 {
2285 /* === Local variables ================================================== */
2286
2287 Int c ; /* a column index */
2288 Int r, row ; /* a row index */
2289 Int *cp ; /* a column pointer */
2290 Int deg ; /* degree of a row or column */
2291 Int *cp_end ; /* a pointer to the end of a column */
2292 Int *new_cp ; /* new column pointer */
2293 Int col_length ; /* length of pruned column */
2294 Int score ; /* current column score */
2295 Int n_col2 ; /* number of non-dense, non-empty columns */
2296 Int n_row2 ; /* number of non-dense, non-empty rows */
2297 Int dense_row_count ; /* remove rows with more entries than this */
2298 Int dense_col_count ; /* remove cols with more entries than this */
2299 Int max_deg ; /* maximum row degree */
2300 Int s ; /* a cset index */
2301 Int ndense_row ; /* number of dense rows */
2302 Int nempty_row ; /* number of empty rows */
2303 Int nnewlyempty_row ; /* number of newly empty rows */
2304 Int ndense_col ; /* number of dense cols (excl "empty" cols) */
2305 Int nempty_col ; /* number of original empty cols */
2306 Int nnewlyempty_col ; /* number of newly empty cols */
2307 Int ne ;
2308
2309 #ifndef NDEBUG
2310 Int debug_count ; /* debug only. */
2311 #endif
2312
2313 /* === Extract knobs ==================================================== */
2314
2315 /* Note: if knobs contains a NaN, this is undefined: */
2316 if (knobs [CCOLAMD_DENSE_ROW] < 0)
2317 {
2318 /* only remove completely dense rows */
2319 dense_row_count = n_col-1 ;
2320 }
2321 else
2322 {
2323 dense_row_count = DENSE_DEGREE (knobs [CCOLAMD_DENSE_ROW], n_col) ;
2324 }
2325 if (knobs [CCOLAMD_DENSE_COL] < 0)
2326 {
2327 /* only remove completely dense columns */
2328 dense_col_count = n_row-1 ;
2329 }
2330 else
2331 {
2332 dense_col_count =
2333 DENSE_DEGREE (knobs [CCOLAMD_DENSE_COL], MIN (n_row, n_col)) ;
2334 }
2335
2336 DEBUG1 (("densecount: "ID" "ID"\n", dense_row_count, dense_col_count)) ;
2337 max_deg = 0 ;
2338
2339 n_col2 = n_col ;
2340 n_row2 = n_row ;
2341
2342 /* Set the head array for bookkeeping of dense and empty columns. */
2343 /* This will be used as hash buckets later. */
2344 for (s = 0 ; s < n_cset ; s++)
2345 {
2346 head [s] = cset_start [s+1] ;
2347 }
2348
2349 ndense_col = 0 ;
2350 nempty_col = 0 ;
2351 nnewlyempty_col = 0 ;
2352 ndense_row = 0 ;
2353 nempty_row = 0 ;
2354 nnewlyempty_row = 0 ;
2355
2356 /* === Kill empty columns =============================================== */
2357
2358 /* Put the empty columns at the end in their natural order, so that LU */
2359 /* factorization can proceed as far as possible. */
2360 for (c = n_col-1 ; c >= 0 ; c--)
2361 {
2362 deg = Col [c].length ;
2363 if (deg == 0)
2364 {
2365 /* this is a empty column, kill and order it last of its cset */
2366 Col [c].shared2.order = --head [CMEMBER (c)] ;
2367 --n_col2 ;
2368 dead_cols [CMEMBER (c)] ++ ;
2369 nempty_col++ ;
2370 KILL_PRINCIPAL_COL (c) ;
2371 }
2372 }
2373 DEBUG1 (("ccolamd: null columns killed: "ID"\n", n_col - n_col2)) ;
2374
2375 /* === Kill dense columns =============================================== */
2376
2377 /* Put the dense columns at the end, in their natural order */
2378 for (c = n_col-1 ; c >= 0 ; c--)
2379 {
2380 /* skip any dead columns */
2381 if (COL_IS_DEAD (c))
2382 {
2383 continue ;
2384 }
2385 deg = Col [c].length ;
2386 if (deg > dense_col_count)
2387 {
2388 /* this is a dense column, kill and order it last of its cset */
2389 Col [c].shared2.order = --head [CMEMBER (c)] ;
2390 --n_col2 ;
2391 dead_cols [CMEMBER (c)] ++ ;
2392 ndense_col++ ;
2393 /* decrement the row degrees */
2394 cp = &A [Col [c].start] ;
2395 cp_end = cp + Col [c].length ;
2396 while (cp < cp_end)
2397 {
2398 Row [*cp++].shared1.degree-- ;
2399 }
2400 KILL_PRINCIPAL_COL (c) ;
2401 }
2402 }
2403 DEBUG1 (("Dense and null columns killed: "ID"\n", n_col - n_col2)) ;
2404
2405 /* === Kill dense and empty rows ======================================== */
2406
2407 /* Note that there can now be empty rows, since dense columns have
2408 * been deleted. These are "newly" empty rows. */
2409
2410 ne = 0 ;
2411 for (r = 0 ; r < n_row ; r++)
2412 {
2413 deg = Row [r].shared1.degree ;
2414 ASSERT (deg >= 0 && deg <= n_col) ;
2415 if (deg > dense_row_count)
2416 {
2417 /* There is at least one dense row. Continue ordering, but */
2418 /* symbolic factorization will be redone after ccolamd is done.*/
2419 ndense_row++ ;
2420 }
2421 if (deg == 0)
2422 {
2423 /* this is a newly empty row, or original empty row */
2424 ne++ ;
2425 }
2426 if (deg > dense_row_count || deg == 0)
2427 {
2428 /* kill a dense or empty row */
2429 KILL_ROW (r) ;
2430 Row [r].thickness = 0 ;
2431 --n_row2 ;
2432 }
2433 else
2434 {
2435 /* keep track of max degree of remaining rows */
2436 max_deg = MAX (max_deg, deg) ;
2437 }
2438 }
2439 nnewlyempty_row = ne - nempty_row ;
2440 DEBUG1 (("ccolamd: Dense and null rows killed: "ID"\n", n_row - n_row2)) ;
2441
2442 /* === Compute initial column scores ==================================== */
2443
2444 /* At this point the row degrees are accurate. They reflect the number */
2445 /* of "live" (non-dense) columns in each row. No empty rows exist. */
2446 /* Some "live" columns may contain only dead rows, however. These are */
2447 /* pruned in the code below. */
2448
2449 /* now find the initial COLMMD score for each column */
2450 for (c = n_col-1 ; c >= 0 ; c--)
2451 {
2452 /* skip dead column */
2453 if (COL_IS_DEAD (c))
2454 {
2455 continue ;
2456 }
2457 score = 0 ;
2458 cp = &A [Col [c].start] ;
2459 new_cp = cp ;
2460 cp_end = cp + Col [c].length ;
2461 while (cp < cp_end)
2462 {
2463 /* get a row */
2464 row = *cp++ ;
2465 /* skip if dead */
2466 if (ROW_IS_DEAD (row))
2467 {
2468 continue ;
2469 }
2470 /* compact the column */
2471 *new_cp++ = row ;
2472 /* add row's external degree */
2473 score += Row [row].shared1.degree - 1 ;
2474 /* guard against integer overflow */
2475 score = MIN (score, n_col) ;
2476 }
2477 /* determine pruned column length */
2478 col_length = (Int) (new_cp - &A [Col [c].start]) ;
2479 if (col_length == 0)
2480 {
2481 /* a newly-made null column (all rows in this col are "dense" */
2482 /* and have already been killed) */
2483 DEBUG1 (("Newly null killed: "ID"\n", c)) ;
2484 Col [c].shared2.order = -- head [CMEMBER (c)] ;
2485 --n_col2 ;
2486 dead_cols [CMEMBER (c)] ++ ;
2487 nnewlyempty_col++ ;
2488 KILL_PRINCIPAL_COL (c) ;
2489 }
2490 else
2491 {
2492 /* set column length and set score */
2493 ASSERT (score >= 0) ;
2494 ASSERT (score <= n_col) ;
2495 Col [c].length = col_length ;
2496 Col [c].shared2.score = score ;
2497 }
2498 }
2499 DEBUG1 (("ccolamd: Dense, null, and newly-null columns killed: "ID"\n",
2500 n_col-n_col2)) ;
2501
2502 /* At this point, all empty rows and columns are dead. All live columns */
2503 /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
2504 /* yet). Rows may contain dead columns, but all live rows contain at */
2505 /* least one live column. */
2506
2507 #ifndef NDEBUG
2508 debug_count = 0 ;
2509 #endif
2510
2511 /* clear the hash buckets */
2512 for (c = 0 ; c <= n_col ; c++)
2513 {
2514 head [c] = EMPTY ;
2515 }
2516
2517 #ifndef NDEBUG
2518 debug_structures (n_row, n_col, Row, Col, A, cmember, cset_start) ;
2519 #endif
2520
2521 /* === Return number of remaining columns, and max row degree =========== */
2522
2523 *p_n_col2 = n_col2 ;
2524 *p_n_row2 = n_row2 ;
2525 *p_max_deg = max_deg ;
2526 *p_ndense_row = ndense_row ;
2527 *p_nempty_row = nempty_row ; /* original empty rows */
2528 *p_nnewlyempty_row = nnewlyempty_row ;
2529 *p_ndense_col = ndense_col ;
2530 *p_nempty_col = nempty_col ; /* original empty cols */
2531 *p_nnewlyempty_col = nnewlyempty_col ;
2532 }
2533
2534
2535 /* ========================================================================== */
2536 /* === find_ordering ======================================================== */
2537 /* ========================================================================== */
2538
2539 /*
2540 * Order the principal columns of the supercolumn form of the matrix
2541 * (no supercolumns on input). Uses a minimum approximate column minimum
2542 * degree ordering method. Not user-callable.
2543 */
2544
find_ordering(Int n_row,Int n_col,Int Alen,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int head[],Int n_col2,Int max_deg,Int pfree,Int cset[],Int cset_start[],Int n_cset,Int cmember[],Int Front_npivcol[],Int Front_nrows[],Int Front_ncols[],Int Front_parent[],Int Front_cols[],Int * p_nfr,Int aggressive,Int InFront[],Int order_for_lu)2545 PRIVATE Int find_ordering /* return the number of garbage collections */
2546 (
2547 /* === Parameters ======================================================= */
2548
2549 Int n_row, /* number of rows of A */
2550 Int n_col, /* number of columns of A */
2551 Int Alen, /* size of A, 2*nnz + n_col or larger */
2552 CColamd_Row Row [ ], /* of size n_row+1 */
2553 CColamd_Col Col [ ], /* of size n_col+1 */
2554 Int A [ ], /* column form and row form of A */
2555 Int head [ ], /* of size n_col+1 */
2556 #ifndef NDEBUG
2557 Int n_col2, /* Remaining columns to order */
2558 #endif
2559 Int max_deg, /* Maximum row degree */
2560 Int pfree, /* index of first free slot (2*nnz on entry) */
2561 Int cset [ ], /* constraint set of A */
2562 Int cset_start [ ], /* pointer to the start of every cset */
2563 #ifndef NDEBUG
2564 Int n_cset, /* number of csets */
2565 #endif
2566 Int cmember [ ], /* col -> cset mapping */
2567 Int Front_npivcol [ ],
2568 Int Front_nrows [ ],
2569 Int Front_ncols [ ],
2570 Int Front_parent [ ],
2571 Int Front_cols [ ],
2572 Int *p_nfr, /* number of fronts */
2573 Int aggressive,
2574 Int InFront [ ],
2575 Int order_for_lu
2576 )
2577 {
2578 /* === Local variables ================================================== */
2579
2580 Int k ; /* current pivot ordering step */
2581 Int pivot_col ; /* current pivot column */
2582 Int *cp ; /* a column pointer */
2583 Int *rp ; /* a row pointer */
2584 Int pivot_row ; /* current pivot row */
2585 Int *new_cp ; /* modified column pointer */
2586 Int *new_rp ; /* modified row pointer */
2587 Int pivot_row_start ; /* pointer to start of pivot row */
2588 Int pivot_row_degree ; /* number of columns in pivot row */
2589 Int pivot_row_length ; /* number of supercolumns in pivot row */
2590 Int pivot_col_score ; /* score of pivot column */
2591 Int needed_memory ; /* free space needed for pivot row */
2592 Int *cp_end ; /* pointer to the end of a column */
2593 Int *rp_end ; /* pointer to the end of a row */
2594 Int row ; /* a row index */
2595 Int col ; /* a column index */
2596 Int max_score ; /* maximum possible score */
2597 Int cur_score ; /* score of current column */
2598 unsigned Int hash ; /* hash value for supernode detection */
2599 Int head_column ; /* head of hash bucket */
2600 Int first_col ; /* first column in hash bucket */
2601 Int tag_mark ; /* marker value for mark array */
2602 Int row_mark ; /* Row [row].shared2.mark */
2603 Int set_difference ; /* set difference size of row with pivot row */
2604 Int min_score ; /* smallest column score */
2605 Int col_thickness ; /* "thickness" (no. of columns in a supercol) */
2606 Int max_mark ; /* maximum value of tag_mark */
2607 Int pivot_col_thickness ; /* number of columns represented by pivot col */
2608 Int prev_col ; /* Used by Dlist operations. */
2609 Int next_col ; /* Used by Dlist operations. */
2610 Int ngarbage ; /* number of garbage collections performed */
2611 Int current_set ; /* consraint set that is being ordered */
2612 Int score ; /* score of a column */
2613 Int colstart ; /* pointer to first column in current cset */
2614 Int colend ; /* pointer to last column in current cset */
2615 Int deadcol ; /* number of dense & null columns in a cset */
2616
2617 #ifndef NDEBUG
2618 Int debug_d ; /* debug loop counter */
2619 Int debug_step = 0 ; /* debug loop counter */
2620 Int cols_thickness = 0 ; /* the thickness of the columns in current */
2621 /* cset degreelist and in pivot row pattern. */
2622 #endif
2623
2624 Int pivot_row_thickness ; /* number of rows represented by pivot row */
2625 Int nfr = 0 ; /* number of fronts */
2626 Int child ;
2627
2628 /* === Initialization and clear mark ==================================== */
2629
2630 max_mark = Int_MAX - n_col ; /* Int_MAX defined in <limits.h> */
2631 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2632 min_score = 0 ;
2633 ngarbage = 0 ;
2634 current_set = -1 ;
2635 deadcol = 0 ;
2636 DEBUG1 (("ccolamd: Ordering, n_col2="ID"\n", n_col2)) ;
2637
2638 for (row = 0 ; row < n_row ; row++)
2639 {
2640 InFront [row] = EMPTY ;
2641 }
2642
2643 /* === Order the columns ================================================ */
2644
2645 for (k = 0 ; k < n_col ; /* 'k' is incremented below */)
2646 {
2647
2648 /* make sure degree list isn't empty */
2649 ASSERT (min_score >= 0) ;
2650 ASSERT (min_score <= n_col) ;
2651 ASSERT (head [min_score] >= EMPTY) ;
2652
2653 #ifndef NDEBUG
2654 for (debug_d = 0 ; debug_d < min_score ; debug_d++)
2655 {
2656 ASSERT (head [debug_d] == EMPTY) ;
2657 }
2658 #endif
2659
2660 /* Initialize the degree list with columns from next non-empty cset */
2661
2662 while ((k+deadcol) == cset_start [current_set+1])
2663 {
2664 current_set++ ;
2665 DEBUG1 (("\n\n\n============ CSET: "ID"\n", current_set)) ;
2666 k += deadcol ; /* jump to start of next cset */
2667 deadcol = 0 ; /* reset dead column count */
2668
2669 ASSERT ((current_set == n_cset) == (k == n_col)) ;
2670
2671 /* return if all columns are ordered. */
2672 if (k == n_col)
2673 {
2674 *p_nfr = nfr ;
2675 return (ngarbage) ;
2676 }
2677
2678 #ifndef NDEBUG
2679 for (col = 0 ; col <= n_col ; col++)
2680 {
2681 ASSERT (head [col] == EMPTY) ;
2682 }
2683 #endif
2684
2685 min_score = n_col ;
2686 colstart = cset_start [current_set] ;
2687 colend = cset_start [current_set+1] ;
2688
2689 while (colstart < colend)
2690 {
2691 col = cset [colstart++] ;
2692
2693 if (COL_IS_DEAD(col))
2694 {
2695 DEBUG1 (("Column "ID" is dead\n", col)) ;
2696 /* count dense and null columns */
2697 if (Col [col].shared2.order != EMPTY)
2698 {
2699 deadcol++ ;
2700 }
2701 continue ;
2702 }
2703
2704 /* only add principal columns in current set to degree lists */
2705 ASSERT (CMEMBER (col) == current_set) ;
2706
2707 score = Col [col].shared2.score ;
2708 DEBUG1 (("Column "ID" is alive, score "ID"\n", col, score)) ;
2709
2710 ASSERT (min_score >= 0) ;
2711 ASSERT (min_score <= n_col) ;
2712 ASSERT (score >= 0) ;
2713 ASSERT (score <= n_col) ;
2714 ASSERT (head [score] >= EMPTY) ;
2715
2716 /* now add this column to dList at proper score location */
2717 next_col = head [score] ;
2718 Col [col].shared3.prev = EMPTY ;
2719 Col [col].shared4.degree_next = next_col ;
2720
2721 /* if there already was a column with the same score, set its */
2722 /* previous pointer to this new column */
2723 if (next_col != EMPTY)
2724 {
2725 Col [next_col].shared3.prev = col ;
2726 }
2727 head [score] = col ;
2728
2729 /* see if this score is less than current min */
2730 min_score = MIN (min_score, score) ;
2731 }
2732
2733 #ifndef NDEBUG
2734 DEBUG1 (("degree lists initialized \n")) ;
2735 debug_deg_lists (n_row, n_col, Row, Col, head, min_score,
2736 ((cset_start [current_set+1]-cset_start [current_set])-deadcol),
2737 max_deg) ;
2738 #endif
2739 }
2740
2741 #ifndef NDEBUG
2742 if (debug_step % 100 == 0)
2743 {
2744 DEBUG2 (("\n... Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
2745 }
2746 else
2747 {
2748 DEBUG3 (("\n------Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
2749 }
2750 debug_step++ ;
2751 DEBUG1 (("start of step k="ID": ", k)) ;
2752 debug_deg_lists (n_row, n_col, Row, Col, head,
2753 min_score, cset_start [current_set+1]-(k+deadcol), max_deg) ;
2754 debug_matrix (n_row, n_col, Row, Col, A) ;
2755 #endif
2756
2757 /* === Select pivot column, and order it ============================ */
2758
2759 while (head [min_score] == EMPTY && min_score < n_col)
2760 {
2761 min_score++ ;
2762 }
2763
2764 pivot_col = head [min_score] ;
2765
2766 ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
2767 next_col = Col [pivot_col].shared4.degree_next ;
2768 head [min_score] = next_col ;
2769 if (next_col != EMPTY)
2770 {
2771 Col [next_col].shared3.prev = EMPTY ;
2772 }
2773
2774 ASSERT (COL_IS_ALIVE (pivot_col)) ;
2775
2776 /* remember score for defrag check */
2777 pivot_col_score = Col [pivot_col].shared2.score ;
2778
2779 /* the pivot column is the kth column in the pivot order */
2780 Col [pivot_col].shared2.order = k ;
2781
2782 /* increment order count by column thickness */
2783 pivot_col_thickness = Col [pivot_col].shared1.thickness ;
2784 k += pivot_col_thickness ;
2785 ASSERT (pivot_col_thickness > 0) ;
2786 DEBUG3 (("Pivot col: "ID" thick "ID"\n", pivot_col,
2787 pivot_col_thickness)) ;
2788
2789 /* === Garbage_collection, if necessary ============================= */
2790
2791 needed_memory = MIN (pivot_col_score, n_col - k) ;
2792 if (pfree + needed_memory >= Alen)
2793 {
2794 pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
2795 ngarbage++ ;
2796 /* after garbage collection we will have enough */
2797 ASSERT (pfree + needed_memory < Alen) ;
2798 /* garbage collection has wiped out Row [ ].shared2.mark array */
2799 tag_mark = clear_mark (0, max_mark, n_row, Row) ;
2800
2801 #ifndef NDEBUG
2802 debug_matrix (n_row, n_col, Row, Col, A) ;
2803 #endif
2804 }
2805
2806 /* === Compute pivot row pattern ==================================== */
2807
2808 /* get starting location for this new merged row */
2809 pivot_row_start = pfree ;
2810
2811 /* initialize new row counts to zero */
2812 pivot_row_degree = 0 ;
2813 pivot_row_thickness = 0 ;
2814
2815 /* tag pivot column as having been visited so it isn't included */
2816 /* in merged pivot row */
2817 Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
2818
2819 /* pivot row is the union of all rows in the pivot column pattern */
2820 cp = &A [Col [pivot_col].start] ;
2821 cp_end = cp + Col [pivot_col].length ;
2822 while (cp < cp_end)
2823 {
2824 /* get a row */
2825 row = *cp++ ;
2826 ASSERT (row >= 0 && row < n_row) ;
2827 DEBUG4 (("Pivcol pattern "ID" "ID"\n", ROW_IS_ALIVE (row), row)) ;
2828 /* skip if row is dead */
2829 if (ROW_IS_ALIVE (row))
2830 {
2831 /* sum the thicknesses of all the rows */
2832 pivot_row_thickness += Row [row].thickness ;
2833
2834 rp = &A [Row [row].start] ;
2835 rp_end = rp + Row [row].length ;
2836 while (rp < rp_end)
2837 {
2838 /* get a column */
2839 col = *rp++ ;
2840 /* add the column, if alive and untagged */
2841 col_thickness = Col [col].shared1.thickness ;
2842 if (col_thickness > 0 && COL_IS_ALIVE (col))
2843 {
2844 /* tag column in pivot row */
2845 Col [col].shared1.thickness = -col_thickness ;
2846 ASSERT (pfree < Alen) ;
2847 /* place column in pivot row */
2848 A [pfree++] = col ;
2849 pivot_row_degree += col_thickness ;
2850 DEBUG4 (("\t\t\tNew live col in pivrow: "ID"\n",col)) ;
2851 }
2852 #ifndef NDEBUG
2853 if (col_thickness < 0 && COL_IS_ALIVE (col))
2854 {
2855 DEBUG4 (("\t\t\tOld live col in pivrow: "ID"\n",col)) ;
2856 }
2857 #endif
2858 }
2859 }
2860 }
2861
2862 /* pivot_row_thickness is the number of rows in frontal matrix */
2863 /* including both pivotal rows and nonpivotal rows */
2864
2865 /* clear tag on pivot column */
2866 Col [pivot_col].shared1.thickness = pivot_col_thickness ;
2867 max_deg = MAX (max_deg, pivot_row_degree) ;
2868
2869 #ifndef NDEBUG
2870 DEBUG3 (("check2\n")) ;
2871 debug_mark (n_row, Row, tag_mark, max_mark) ;
2872 #endif
2873
2874 /* === Kill all rows used to construct pivot row ==================== */
2875
2876 /* also kill pivot row, temporarily */
2877 cp = &A [Col [pivot_col].start] ;
2878 cp_end = cp + Col [pivot_col].length ;
2879 while (cp < cp_end)
2880 {
2881 /* may be killing an already dead row */
2882 row = *cp++ ;
2883 DEBUG3 (("Kill row in pivot col: "ID"\n", row)) ;
2884 ASSERT (row >= 0 && row < n_row) ;
2885 if (ROW_IS_ALIVE (row))
2886 {
2887 if (Row [row].front != EMPTY)
2888 {
2889 /* This row represents a frontal matrix. */
2890 /* Row [row].front is a child of current front */
2891 child = Row [row].front ;
2892 Front_parent [child] = nfr ;
2893 DEBUG1 (("Front "ID" => front "ID", normal\n", child, nfr));
2894 }
2895 else
2896 {
2897 /* This is an original row. Keep track of which front
2898 * is its parent in the row-merge tree. */
2899 InFront [row] = nfr ;
2900 DEBUG1 (("Row "ID" => front "ID", normal\n", row, nfr)) ;
2901 }
2902 }
2903
2904 KILL_ROW (row) ;
2905 Row [row].thickness = 0 ;
2906 }
2907
2908 /* === Select a row index to use as the new pivot row =============== */
2909
2910 pivot_row_length = pfree - pivot_row_start ;
2911 if (pivot_row_length > 0)
2912 {
2913 /* pick the "pivot" row arbitrarily (first row in col) */
2914 pivot_row = A [Col [pivot_col].start] ;
2915 DEBUG3 (("Pivotal row is "ID"\n", pivot_row)) ;
2916 }
2917 else
2918 {
2919 /* there is no pivot row, since it is of zero length */
2920 pivot_row = EMPTY ;
2921 ASSERT (pivot_row_length == 0) ;
2922 }
2923 ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
2924
2925 /* === Approximate degree computation =============================== */
2926
2927 /* Here begins the computation of the approximate degree. The column */
2928 /* score is the sum of the pivot row "length", plus the size of the */
2929 /* set differences of each row in the column minus the pattern of the */
2930 /* pivot row itself. The column ("thickness") itself is also */
2931 /* excluded from the column score (we thus use an approximate */
2932 /* external degree). */
2933
2934 /* The time taken by the following code (compute set differences, and */
2935 /* add them up) is proportional to the size of the data structure */
2936 /* being scanned - that is, the sum of the sizes of each column in */
2937 /* the pivot row. Thus, the amortized time to compute a column score */
2938 /* is proportional to the size of that column (where size, in this */
2939 /* context, is the column "length", or the number of row indices */
2940 /* in that column). The number of row indices in a column is */
2941 /* monotonically non-decreasing, from the length of the original */
2942 /* column on input to colamd. */
2943
2944 /* === Compute set differences ====================================== */
2945
2946 DEBUG3 (("** Computing set differences phase. **\n")) ;
2947
2948 /* pivot row is currently dead - it will be revived later. */
2949
2950 DEBUG3 (("Pivot row: ")) ;
2951 /* for each column in pivot row */
2952 rp = &A [pivot_row_start] ;
2953 rp_end = rp + pivot_row_length ;
2954 while (rp < rp_end)
2955 {
2956 col = *rp++ ;
2957 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
2958 DEBUG3 (("Col: "ID"\n", col)) ;
2959
2960 /* clear tags used to construct pivot row pattern */
2961 col_thickness = -Col [col].shared1.thickness ;
2962 ASSERT (col_thickness > 0) ;
2963 Col [col].shared1.thickness = col_thickness ;
2964
2965 /* === Remove column from degree list =========================== */
2966
2967 /* only columns in current_set will be in degree list */
2968 if (CMEMBER (col) == current_set)
2969 {
2970 #ifndef NDEBUG
2971 cols_thickness += col_thickness ;
2972 #endif
2973 cur_score = Col [col].shared2.score ;
2974 prev_col = Col [col].shared3.prev ;
2975 next_col = Col [col].shared4.degree_next ;
2976 DEBUG3 ((" cur_score "ID" prev_col "ID" next_col "ID"\n",
2977 cur_score, prev_col, next_col)) ;
2978 ASSERT (cur_score >= 0) ;
2979 ASSERT (cur_score <= n_col) ;
2980 ASSERT (cur_score >= EMPTY) ;
2981 if (prev_col == EMPTY)
2982 {
2983 head [cur_score] = next_col ;
2984 }
2985 else
2986 {
2987 Col [prev_col].shared4.degree_next = next_col ;
2988 }
2989 if (next_col != EMPTY)
2990 {
2991 Col [next_col].shared3.prev = prev_col ;
2992 }
2993 }
2994
2995 /* === Scan the column ========================================== */
2996
2997 cp = &A [Col [col].start] ;
2998 cp_end = cp + Col [col].length ;
2999 while (cp < cp_end)
3000 {
3001 /* get a row */
3002 row = *cp++ ;
3003 row_mark = Row [row].shared2.mark ;
3004 /* skip if dead */
3005 if (ROW_IS_MARKED_DEAD (row_mark))
3006 {
3007 continue ;
3008 }
3009 ASSERT (row != pivot_row) ;
3010 set_difference = row_mark - tag_mark ;
3011 /* check if the row has been seen yet */
3012 if (set_difference < 0)
3013 {
3014 ASSERT (Row [row].shared1.degree <= max_deg) ;
3015 set_difference = Row [row].shared1.degree ;
3016 }
3017 /* subtract column thickness from this row's set difference */
3018 set_difference -= col_thickness ;
3019 ASSERT (set_difference >= 0) ;
3020 /* absorb this row if the set difference becomes zero */
3021 if (set_difference == 0 && aggressive)
3022 {
3023 DEBUG3 (("aggressive absorption. Row: "ID"\n", row)) ;
3024
3025 if (Row [row].front != EMPTY)
3026 {
3027 /* Row [row].front is a child of current front. */
3028 child = Row [row].front ;
3029 Front_parent [child] = nfr ;
3030 DEBUG1 (("Front "ID" => front "ID", aggressive\n",
3031 child, nfr)) ;
3032 }
3033 else
3034 {
3035 /* this is an original row. Keep track of which front
3036 * assembles it, for the row-merge tree */
3037 InFront [row] = nfr ;
3038 DEBUG1 (("Row "ID" => front "ID", aggressive\n",
3039 row, nfr)) ;
3040 }
3041
3042 KILL_ROW (row) ;
3043
3044 /* sum the thicknesses of all the rows */
3045 pivot_row_thickness += Row [row].thickness ;
3046 Row [row].thickness = 0 ;
3047 }
3048 else
3049 {
3050 /* save the new mark */
3051 Row [row].shared2.mark = set_difference + tag_mark ;
3052 }
3053 }
3054 }
3055
3056 #ifndef NDEBUG
3057 debug_deg_lists (n_row, n_col, Row, Col, head, min_score,
3058 cset_start [current_set+1]-(k+deadcol)-(cols_thickness),
3059 max_deg) ;
3060 cols_thickness = 0 ;
3061 #endif
3062
3063 /* === Add up set differences for each column ======================= */
3064
3065 DEBUG3 (("** Adding set differences phase. **\n")) ;
3066
3067 /* for each column in pivot row */
3068 rp = &A [pivot_row_start] ;
3069 rp_end = rp + pivot_row_length ;
3070 while (rp < rp_end)
3071 {
3072 /* get a column */
3073 col = *rp++ ;
3074 ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
3075 hash = 0 ;
3076 cur_score = 0 ;
3077 cp = &A [Col [col].start] ;
3078 /* compact the column */
3079 new_cp = cp ;
3080 cp_end = cp + Col [col].length ;
3081
3082 DEBUG4 (("Adding set diffs for Col: "ID".\n", col)) ;
3083
3084 while (cp < cp_end)
3085 {
3086 /* get a row */
3087 row = *cp++ ;
3088 ASSERT (row >= 0 && row < n_row) ;
3089 row_mark = Row [row].shared2.mark ;
3090 /* skip if dead */
3091 if (ROW_IS_MARKED_DEAD (row_mark))
3092 {
3093 DEBUG4 ((" Row "ID", dead\n", row)) ;
3094 continue ;
3095 }
3096 DEBUG4 ((" Row "ID", set diff "ID"\n", row, row_mark-tag_mark));
3097 ASSERT (row_mark >= tag_mark) ;
3098 /* compact the column */
3099 *new_cp++ = row ;
3100 /* compute hash function */
3101 hash += row ;
3102 /* add set difference */
3103 cur_score += row_mark - tag_mark ;
3104 /* integer overflow... */
3105 cur_score = MIN (cur_score, n_col) ;
3106 }
3107
3108 /* recompute the column's length */
3109 Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
3110
3111 /* === Further mass elimination ================================= */
3112
3113 if (Col [col].length == 0 && CMEMBER (col) == current_set)
3114 {
3115 DEBUG4 (("further mass elimination. Col: "ID"\n", col)) ;
3116 /* nothing left but the pivot row in this column */
3117 KILL_PRINCIPAL_COL (col) ;
3118 pivot_row_degree -= Col [col].shared1.thickness ;
3119 ASSERT (pivot_row_degree >= 0) ;
3120 /* order it */
3121 Col [col].shared2.order = k ;
3122 /* increment order count by column thickness */
3123 k += Col [col].shared1.thickness ;
3124 pivot_col_thickness += Col [col].shared1.thickness ;
3125 /* add to column list of front */
3126 #ifndef NDEBUG
3127 DEBUG1 (("Mass")) ;
3128 dump_super (col, Col, n_col) ;
3129 #endif
3130 Col [Col [col].lastcol].nextcol = Front_cols [nfr] ;
3131 Front_cols [nfr] = col ;
3132 }
3133 else
3134 {
3135 /* === Prepare for supercolumn detection ==================== */
3136
3137 DEBUG4 (("Preparing supercol detection for Col: "ID".\n", col));
3138
3139 /* save score so far */
3140 Col [col].shared2.score = cur_score ;
3141
3142 /* add column to hash table, for supercolumn detection */
3143 hash %= n_col + 1 ;
3144
3145 DEBUG4 ((" Hash = "ID", n_col = "ID".\n", hash, n_col)) ;
3146 ASSERT (((Int) hash) <= n_col) ;
3147
3148 head_column = head [hash] ;
3149 if (head_column > EMPTY)
3150 {
3151 /* degree list "hash" is non-empty, use prev (shared3) of */
3152 /* first column in degree list as head of hash bucket */
3153 first_col = Col [head_column].shared3.headhash ;
3154 Col [head_column].shared3.headhash = col ;
3155 }
3156 else
3157 {
3158 /* degree list "hash" is empty, use head as hash bucket */
3159 first_col = - (head_column + 2) ;
3160 head [hash] = - (col + 2) ;
3161 }
3162 Col [col].shared4.hash_next = first_col ;
3163
3164 /* save hash function in Col [col].shared3.hash */
3165 Col [col].shared3.hash = (Int) hash ;
3166 ASSERT (COL_IS_ALIVE (col)) ;
3167 }
3168 }
3169
3170 /* The approximate external column degree is now computed. */
3171
3172 /* === Supercolumn detection ======================================== */
3173
3174 DEBUG3 (("** Supercolumn detection phase. **\n")) ;
3175
3176 detect_super_cols (
3177 #ifndef NDEBUG
3178 n_col, Row,
3179 #endif
3180 Col, A, head, pivot_row_start, pivot_row_length, cmember) ;
3181
3182 /* === Kill the pivotal column ====================================== */
3183
3184 DEBUG1 ((" KILLING column detect supercols "ID" \n", pivot_col)) ;
3185 KILL_PRINCIPAL_COL (pivot_col) ;
3186
3187 /* add columns to column list of front */
3188 #ifndef NDEBUG
3189 DEBUG1 (("Pivot")) ;
3190 dump_super (pivot_col, Col, n_col) ;
3191 #endif
3192 Col [Col [pivot_col].lastcol].nextcol = Front_cols [nfr] ;
3193 Front_cols [nfr] = pivot_col ;
3194
3195 /* === Clear mark =================================================== */
3196
3197 tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
3198
3199 #ifndef NDEBUG
3200 DEBUG3 (("check3\n")) ;
3201 debug_mark (n_row, Row, tag_mark, max_mark) ;
3202 #endif
3203
3204 /* === Finalize the new pivot row, and column scores ================ */
3205
3206 DEBUG3 (("** Finalize scores phase. **\n")) ;
3207
3208 /* for each column in pivot row */
3209 rp = &A [pivot_row_start] ;
3210 /* compact the pivot row */
3211 new_rp = rp ;
3212 rp_end = rp + pivot_row_length ;
3213 while (rp < rp_end)
3214 {
3215 col = *rp++ ;
3216 /* skip dead columns */
3217 if (COL_IS_DEAD (col))
3218 {
3219 continue ;
3220 }
3221 *new_rp++ = col ;
3222 /* add new pivot row to column */
3223 A [Col [col].start + (Col [col].length++)] = pivot_row ;
3224
3225 /* retrieve score so far and add on pivot row's degree. */
3226 /* (we wait until here for this in case the pivot */
3227 /* row's degree was reduced due to mass elimination). */
3228 cur_score = Col [col].shared2.score + pivot_row_degree ;
3229
3230 /* calculate the max possible score as the number of */
3231 /* external columns minus the 'k' value minus the */
3232 /* columns thickness */
3233 max_score = n_col - k - Col [col].shared1.thickness ;
3234
3235 /* make the score the external degree of the union-of-rows */
3236 cur_score -= Col [col].shared1.thickness ;
3237
3238 /* make sure score is less or equal than the max score */
3239 cur_score = MIN (cur_score, max_score) ;
3240 ASSERT (cur_score >= 0) ;
3241
3242 /* store updated score */
3243 Col [col].shared2.score = cur_score ;
3244
3245 /* === Place column back in degree list ========================= */
3246
3247 if (CMEMBER (col) == current_set)
3248 {
3249 ASSERT (min_score >= 0) ;
3250 ASSERT (min_score <= n_col) ;
3251 ASSERT (cur_score >= 0) ;
3252 ASSERT (cur_score <= n_col) ;
3253 ASSERT (head [cur_score] >= EMPTY) ;
3254 next_col = head [cur_score] ;
3255 Col [col].shared4.degree_next = next_col ;
3256 Col [col].shared3.prev = EMPTY ;
3257 if (next_col != EMPTY)
3258 {
3259 Col [next_col].shared3.prev = col ;
3260 }
3261 head [cur_score] = col ;
3262 /* see if this score is less than current min */
3263 min_score = MIN (min_score, cur_score) ;
3264 }
3265 else
3266 {
3267 Col [col].shared4.degree_next = EMPTY ;
3268 Col [col].shared3.prev = EMPTY ;
3269 }
3270 }
3271
3272 #ifndef NDEBUG
3273 debug_deg_lists (n_row, n_col, Row, Col, head,
3274 min_score, cset_start [current_set+1]-(k+deadcol), max_deg) ;
3275 #endif
3276
3277 /* frontal matrix can have more pivot cols than pivot rows for */
3278 /* singular matrices. */
3279
3280 /* number of candidate pivot columns */
3281 Front_npivcol [nfr] = pivot_col_thickness ;
3282
3283 /* all rows (not just size of contrib. block) */
3284 Front_nrows [nfr] = pivot_row_thickness ;
3285
3286 /* all cols */
3287 Front_ncols [nfr] = pivot_col_thickness + pivot_row_degree ;
3288
3289 Front_parent [nfr] = EMPTY ;
3290
3291 pivot_row_thickness -= pivot_col_thickness ;
3292 DEBUG1 (("Front "ID" Pivot_row_thickness after pivot cols elim: "ID"\n",
3293 nfr, pivot_row_thickness)) ;
3294 pivot_row_thickness = MAX (0, pivot_row_thickness) ;
3295
3296 /* === Resurrect the new pivot row ================================== */
3297
3298 if ((pivot_row_degree > 0 && pivot_row_thickness > 0 && (order_for_lu))
3299 || (pivot_row_degree > 0 && (!order_for_lu)))
3300 {
3301 /* update pivot row length to reflect any cols that were killed */
3302 /* during super-col detection and mass elimination */
3303 Row [pivot_row].start = pivot_row_start ;
3304 Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
3305 Row [pivot_row].shared1.degree = pivot_row_degree ;
3306 Row [pivot_row].shared2.mark = 0 ;
3307 Row [pivot_row].thickness = pivot_row_thickness ;
3308 Row [pivot_row].front = nfr ;
3309 /* pivot row is no longer dead */
3310 DEBUG1 (("Resurrect Pivot_row "ID" deg: "ID"\n",
3311 pivot_row, pivot_row_degree)) ;
3312 }
3313
3314 #ifndef NDEBUG
3315 DEBUG1 (("Front "ID" : "ID" "ID" "ID" ", nfr,
3316 Front_npivcol [nfr], Front_nrows [nfr], Front_ncols [nfr])) ;
3317 DEBUG1 ((" cols:[ ")) ;
3318 debug_d = 0 ;
3319 for (col = Front_cols [nfr] ; col != EMPTY ; col = Col [col].nextcol)
3320 {
3321 DEBUG1 ((" "ID, col)) ;
3322 ASSERT (col >= 0 && col < n_col) ;
3323 ASSERT (COL_IS_DEAD (col)) ;
3324 debug_d++ ;
3325 ASSERT (debug_d <= pivot_col_thickness) ;
3326 }
3327 ASSERT (debug_d == pivot_col_thickness) ;
3328 DEBUG1 ((" ]\n ")) ;
3329 #endif
3330 nfr++ ; /* one more front */
3331 }
3332
3333 /* === All principal columns have now been ordered ====================== */
3334
3335 *p_nfr = nfr ;
3336 return (ngarbage) ;
3337 }
3338
3339
3340 /* ========================================================================== */
3341 /* === detect_super_cols ==================================================== */
3342 /* ========================================================================== */
3343
3344 /*
3345 * Detects supercolumns by finding matches between columns in the hash buckets.
3346 * Check amongst columns in the set A [row_start ... row_start + row_length-1].
3347 * The columns under consideration are currently *not* in the degree lists,
3348 * and have already been placed in the hash buckets.
3349 *
3350 * The hash bucket for columns whose hash function is equal to h is stored
3351 * as follows:
3352 *
3353 * if head [h] is >= 0, then head [h] contains a degree list, so:
3354 *
3355 * head [h] is the first column in degree bucket h.
3356 * Col [head [h]].headhash gives the first column in hash bucket h.
3357 *
3358 * otherwise, the degree list is empty, and:
3359 *
3360 * -(head [h] + 2) is the first column in hash bucket h.
3361 *
3362 * For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
3363 * column" pointer. Col [c].shared3.hash is used instead as the hash number
3364 * for that column. The value of Col [c].shared4.hash_next is the next column
3365 * in the same hash bucket.
3366 *
3367 * Assuming no, or "few" hash collisions, the time taken by this routine is
3368 * linear in the sum of the sizes (lengths) of each column whose score has
3369 * just been computed in the approximate degree computation.
3370 * Not user-callable.
3371 */
3372
detect_super_cols(Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int head[],Int row_start,Int row_length,Int cmember[])3373 PRIVATE void detect_super_cols
3374 (
3375 /* === Parameters ======================================================= */
3376
3377 #ifndef NDEBUG
3378 /* these two parameters are only needed when debugging is enabled: */
3379 Int n_col, /* number of columns of A */
3380 CColamd_Row Row [ ], /* of size n_row+1 */
3381 #endif
3382
3383 CColamd_Col Col [ ], /* of size n_col+1 */
3384 Int A [ ], /* row indices of A */
3385 Int head [ ], /* head of degree lists and hash buckets */
3386 Int row_start, /* pointer to set of columns to check */
3387 Int row_length, /* number of columns to check */
3388 Int cmember [ ] /* col -> cset mapping */
3389 )
3390 {
3391 /* === Local variables ================================================== */
3392
3393 Int hash ; /* hash value for a column */
3394 Int *rp ; /* pointer to a row */
3395 Int c ; /* a column index */
3396 Int super_c ; /* column index of the column to absorb into */
3397 Int *cp1 ; /* column pointer for column super_c */
3398 Int *cp2 ; /* column pointer for column c */
3399 Int length ; /* length of column super_c */
3400 Int prev_c ; /* column preceding c in hash bucket */
3401 Int i ; /* loop counter */
3402 Int *rp_end ; /* pointer to the end of the row */
3403 Int col ; /* a column index in the row to check */
3404 Int head_column ; /* first column in hash bucket or degree list */
3405 Int first_col ; /* first column in hash bucket */
3406
3407 /* === Consider each column in the row ================================== */
3408
3409 rp = &A [row_start] ;
3410 rp_end = rp + row_length ;
3411 while (rp < rp_end)
3412 {
3413 col = *rp++ ;
3414 if (COL_IS_DEAD (col))
3415 {
3416 continue ;
3417 }
3418
3419 /* get hash number for this column */
3420 hash = Col [col].shared3.hash ;
3421 ASSERT (hash <= n_col) ;
3422
3423 /* === Get the first column in this hash bucket ===================== */
3424
3425 head_column = head [hash] ;
3426 if (head_column > EMPTY)
3427 {
3428 first_col = Col [head_column].shared3.headhash ;
3429 }
3430 else
3431 {
3432 first_col = - (head_column + 2) ;
3433 }
3434
3435 /* === Consider each column in the hash bucket ====================== */
3436
3437 for (super_c = first_col ; super_c != EMPTY ;
3438 super_c = Col [super_c].shared4.hash_next)
3439 {
3440 ASSERT (COL_IS_ALIVE (super_c)) ;
3441 ASSERT (Col [super_c].shared3.hash == hash) ;
3442 length = Col [super_c].length ;
3443
3444 /* prev_c is the column preceding column c in the hash bucket */
3445 prev_c = super_c ;
3446
3447 /* === Compare super_c with all columns after it ================ */
3448
3449 for (c = Col [super_c].shared4.hash_next ;
3450 c != EMPTY ; c = Col [c].shared4.hash_next)
3451 {
3452 ASSERT (c != super_c) ;
3453 ASSERT (COL_IS_ALIVE (c)) ;
3454 ASSERT (Col [c].shared3.hash == hash) ;
3455
3456 /* not identical if lengths or scores are different, */
3457 /* or if in different constraint sets */
3458 if (Col [c].length != length ||
3459 Col [c].shared2.score != Col [super_c].shared2.score
3460 || CMEMBER (c) != CMEMBER (super_c))
3461 {
3462 prev_c = c ;
3463 continue ;
3464 }
3465
3466 /* compare the two columns */
3467 cp1 = &A [Col [super_c].start] ;
3468 cp2 = &A [Col [c].start] ;
3469
3470 for (i = 0 ; i < length ; i++)
3471 {
3472 /* the columns are "clean" (no dead rows) */
3473 ASSERT (ROW_IS_ALIVE (*cp1)) ;
3474 ASSERT (ROW_IS_ALIVE (*cp2)) ;
3475 /* row indices will same order for both supercols, */
3476 /* no gather scatter nessasary */
3477 if (*cp1++ != *cp2++)
3478 {
3479 break ;
3480 }
3481 }
3482
3483 /* the two columns are different if the for-loop "broke" */
3484 /* super columns should belong to the same constraint set */
3485 if (i != length)
3486 {
3487 prev_c = c ;
3488 continue ;
3489 }
3490
3491 /* === Got it! two columns are identical =================== */
3492
3493 ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
3494
3495 Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
3496 Col [c].shared1.parent = super_c ;
3497 KILL_NON_PRINCIPAL_COL (c) ;
3498 /* order c later, in order_children() */
3499 Col [c].shared2.order = EMPTY ;
3500 /* remove c from hash bucket */
3501 Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
3502
3503 /* add c to end of list of super_c */
3504 ASSERT (Col [super_c].lastcol >= 0) ;
3505 ASSERT (Col [super_c].lastcol < n_col) ;
3506 Col [Col [super_c].lastcol].nextcol = c ;
3507 Col [super_c].lastcol = Col [c].lastcol ;
3508 #ifndef NDEBUG
3509 /* dump the supercolumn */
3510 DEBUG1 (("Super")) ;
3511 dump_super (super_c, Col, n_col) ;
3512 #endif
3513 }
3514 }
3515
3516 /* === Empty this hash bucket ======================================= */
3517
3518 if (head_column > EMPTY)
3519 {
3520 /* corresponding degree list "hash" is not empty */
3521 Col [head_column].shared3.headhash = EMPTY ;
3522 }
3523 else
3524 {
3525 /* corresponding degree list "hash" is empty */
3526 head [hash] = EMPTY ;
3527 }
3528 }
3529 }
3530
3531
3532 /* ========================================================================== */
3533 /* === garbage_collection =================================================== */
3534 /* ========================================================================== */
3535
3536 /*
3537 * Defragments and compacts columns and rows in the workspace A. Used when
3538 * all avaliable memory has been used while performing row merging. Returns
3539 * the index of the first free position in A, after garbage collection. The
3540 * time taken by this routine is linear is the size of the array A, which is
3541 * itself linear in the number of nonzeros in the input matrix.
3542 * Not user-callable.
3543 */
3544
garbage_collection(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int * pfree)3545 PRIVATE Int garbage_collection /* returns the new value of pfree */
3546 (
3547 /* === Parameters ======================================================= */
3548
3549 Int n_row, /* number of rows */
3550 Int n_col, /* number of columns */
3551 CColamd_Row Row [ ], /* row info */
3552 CColamd_Col Col [ ], /* column info */
3553 Int A [ ], /* A [0 ... Alen-1] holds the matrix */
3554 Int *pfree /* &A [0] ... pfree is in use */
3555 )
3556 {
3557 /* === Local variables ================================================== */
3558
3559 Int *psrc ; /* source pointer */
3560 Int *pdest ; /* destination pointer */
3561 Int j ; /* counter */
3562 Int r ; /* a row index */
3563 Int c ; /* a column index */
3564 Int length ; /* length of a row or column */
3565
3566 #ifndef NDEBUG
3567 Int debug_rows ;
3568 DEBUG2 (("Defrag..\n")) ;
3569 for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
3570 debug_rows = 0 ;
3571 #endif
3572
3573 /* === Defragment the columns =========================================== */
3574
3575 pdest = &A[0] ;
3576 for (c = 0 ; c < n_col ; c++)
3577 {
3578 if (COL_IS_ALIVE (c))
3579 {
3580 psrc = &A [Col [c].start] ;
3581
3582 /* move and compact the column */
3583 ASSERT (pdest <= psrc) ;
3584 Col [c].start = (Int) (pdest - &A [0]) ;
3585 length = Col [c].length ;
3586 for (j = 0 ; j < length ; j++)
3587 {
3588 r = *psrc++ ;
3589 if (ROW_IS_ALIVE (r))
3590 {
3591 *pdest++ = r ;
3592 }
3593 }
3594 Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
3595 }
3596 }
3597
3598 /* === Prepare to defragment the rows =================================== */
3599
3600 for (r = 0 ; r < n_row ; r++)
3601 {
3602 if (ROW_IS_DEAD (r) || (Row [r].length == 0))
3603 {
3604 /* This row is already dead, or is of zero length. Cannot compact
3605 * a row of zero length, so kill it. NOTE: in the current version,
3606 * there are no zero-length live rows. Kill the row (for the first
3607 * time, or again) just to be safe. */
3608 KILL_ROW (r) ;
3609 }
3610 else
3611 {
3612 /* save first column index in Row [r].shared2.first_column */
3613 psrc = &A [Row [r].start] ;
3614 Row [r].shared2.first_column = *psrc ;
3615 ASSERT (ROW_IS_ALIVE (r)) ;
3616 /* flag the start of the row with the one's complement of row */
3617 *psrc = ONES_COMPLEMENT (r) ;
3618 #ifndef NDEBUG
3619 debug_rows++ ;
3620 #endif
3621 }
3622 }
3623
3624 /* === Defragment the rows ============================================== */
3625
3626 psrc = pdest ;
3627 while (psrc < pfree)
3628 {
3629 /* find a negative number ... the start of a row */
3630 if (*psrc++ < 0)
3631 {
3632 psrc-- ;
3633 /* get the row index */
3634 r = ONES_COMPLEMENT (*psrc) ;
3635 ASSERT (r >= 0 && r < n_row) ;
3636 /* restore first column index */
3637 *psrc = Row [r].shared2.first_column ;
3638 ASSERT (ROW_IS_ALIVE (r)) ;
3639
3640 /* move and compact the row */
3641 ASSERT (pdest <= psrc) ;
3642 Row [r].start = (Int) (pdest - &A [0]) ;
3643 length = Row [r].length ;
3644 for (j = 0 ; j < length ; j++)
3645 {
3646 c = *psrc++ ;
3647 if (COL_IS_ALIVE (c))
3648 {
3649 *pdest++ = c ;
3650 }
3651 }
3652 Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
3653 #ifndef NDEBUG
3654 debug_rows-- ;
3655 #endif
3656 }
3657 }
3658
3659 /* ensure we found all the rows */
3660 ASSERT (debug_rows == 0) ;
3661
3662 /* === Return the new value of pfree ==================================== */
3663
3664 return ((Int) (pdest - &A [0])) ;
3665 }
3666
3667
3668 /* ========================================================================== */
3669 /* === clear_mark =========================================================== */
3670 /* ========================================================================== */
3671
3672 /*
3673 * Clears the Row [ ].shared2.mark array, and returns the new tag_mark.
3674 * Return value is the new tag_mark. Not user-callable.
3675 */
3676
clear_mark(Int tag_mark,Int max_mark,Int n_row,CColamd_Row Row[])3677 PRIVATE Int clear_mark /* return the new value for tag_mark */
3678 (
3679 /* === Parameters ======================================================= */
3680
3681 Int tag_mark, /* new value of tag_mark */
3682 Int max_mark, /* max allowed value of tag_mark */
3683
3684 Int n_row, /* number of rows in A */
3685 CColamd_Row Row [ ] /* Row [0 ... n_row-1].shared2.mark is set to zero */
3686 )
3687 {
3688 /* === Local variables ================================================== */
3689
3690 Int r ;
3691
3692 if (tag_mark <= 0 || tag_mark >= max_mark)
3693 {
3694 for (r = 0 ; r < n_row ; r++)
3695 {
3696 if (ROW_IS_ALIVE (r))
3697 {
3698 Row [r].shared2.mark = 0 ;
3699 }
3700 }
3701 tag_mark = 1 ;
3702 }
3703
3704 return (tag_mark) ;
3705 }
3706
3707
3708 /* ========================================================================== */
3709 /* === print_report ========================================================= */
3710 /* ========================================================================== */
3711
3712 /* No printing occurs if NPRINT is defined at compile time. */
3713
print_report(char * method,Int stats[CCOLAMD_STATS])3714 PRIVATE void print_report
3715 (
3716 char *method,
3717 Int stats [CCOLAMD_STATS]
3718 )
3719 {
3720
3721 Int i1, i2, i3 ;
3722
3723 SUITESPARSE_PRINTF (("\n%s version %d.%d, %s: ", method,
3724 CCOLAMD_MAIN_VERSION, CCOLAMD_SUB_VERSION, CCOLAMD_DATE)) ;
3725
3726 if (!stats)
3727 {
3728 SUITESPARSE_PRINTF (("No statistics available.\n")) ;
3729 return ;
3730 }
3731
3732 i1 = stats [CCOLAMD_INFO1] ;
3733 i2 = stats [CCOLAMD_INFO2] ;
3734 i3 = stats [CCOLAMD_INFO3] ;
3735
3736 if (stats [CCOLAMD_STATUS] >= 0)
3737 {
3738 SUITESPARSE_PRINTF(("OK. ")) ;
3739 }
3740 else
3741 {
3742 SUITESPARSE_PRINTF(("ERROR. ")) ;
3743 }
3744
3745 switch (stats [CCOLAMD_STATUS])
3746 {
3747
3748 case CCOLAMD_OK_BUT_JUMBLED:
3749
3750 SUITESPARSE_PRINTF((
3751 "Matrix has unsorted or duplicate row indices.\n")) ;
3752
3753 SUITESPARSE_PRINTF((
3754 "%s: duplicate or out-of-order row indices: "ID"\n",
3755 method, i3)) ;
3756
3757 SUITESPARSE_PRINTF((
3758 "%s: last seen duplicate or out-of-order row: "ID"\n",
3759 method, INDEX (i2))) ;
3760
3761 SUITESPARSE_PRINTF((
3762 "%s: last seen in column: "ID"",
3763 method, INDEX (i1))) ;
3764
3765 /* no break - fall through to next case instead */
3766
3767 case CCOLAMD_OK:
3768
3769 SUITESPARSE_PRINTF(("\n")) ;
3770
3771 SUITESPARSE_PRINTF((
3772 "%s: number of dense or empty rows ignored: "ID"\n",
3773 method, stats [CCOLAMD_DENSE_ROW])) ;
3774
3775 SUITESPARSE_PRINTF((
3776 "%s: number of dense or empty columns ignored: "ID"\n",
3777 method, stats [CCOLAMD_DENSE_COL])) ;
3778
3779 SUITESPARSE_PRINTF((
3780 "%s: number of garbage collections performed: "ID"\n",
3781 method, stats [CCOLAMD_DEFRAG_COUNT])) ;
3782 break ;
3783
3784 case CCOLAMD_ERROR_A_not_present:
3785
3786 SUITESPARSE_PRINTF((
3787 "Array A (row indices of matrix) not present.\n")) ;
3788 break ;
3789
3790 case CCOLAMD_ERROR_p_not_present:
3791
3792 SUITESPARSE_PRINTF((
3793 "Array p (column pointers for matrix) not present.\n")) ;
3794 break ;
3795
3796 case CCOLAMD_ERROR_nrow_negative:
3797
3798 SUITESPARSE_PRINTF(("Invalid number of rows ("ID").\n", i1)) ;
3799 break ;
3800
3801 case CCOLAMD_ERROR_ncol_negative:
3802
3803 SUITESPARSE_PRINTF(("Invalid number of columns ("ID").\n", i1)) ;
3804 break ;
3805
3806 case CCOLAMD_ERROR_nnz_negative:
3807
3808 SUITESPARSE_PRINTF((
3809 "Invalid number of nonzero entries ("ID").\n", i1)) ;
3810 break ;
3811
3812 case CCOLAMD_ERROR_p0_nonzero:
3813
3814 SUITESPARSE_PRINTF((
3815 "Invalid column pointer, p [0] = "ID", must be 0.\n", i1)) ;
3816 break ;
3817
3818 case CCOLAMD_ERROR_A_too_small:
3819
3820 SUITESPARSE_PRINTF(("Array A too small.\n")) ;
3821 SUITESPARSE_PRINTF((
3822 " Need Alen >= "ID", but given only Alen = "ID".\n",
3823 i1, i2)) ;
3824 break ;
3825
3826 case CCOLAMD_ERROR_col_length_negative:
3827
3828 SUITESPARSE_PRINTF((
3829 "Column "ID" has a negative number of entries ("ID").\n",
3830 INDEX (i1), i2)) ;
3831 break ;
3832
3833 case CCOLAMD_ERROR_row_index_out_of_bounds:
3834
3835 SUITESPARSE_PRINTF((
3836 "Row index (row "ID") out of bounds ("ID" to "ID") in"
3837 "column "ID".\n", INDEX (i2), INDEX (0), INDEX (i3-1),
3838 INDEX (i1))) ;
3839 break ;
3840
3841 case CCOLAMD_ERROR_out_of_memory:
3842
3843 SUITESPARSE_PRINTF(("Out of memory.\n")) ;
3844 break ;
3845
3846 case CCOLAMD_ERROR_invalid_cmember:
3847
3848 SUITESPARSE_PRINTF(("cmember invalid\n")) ;
3849 break ;
3850 }
3851 }
3852
3853
3854 /* ========================================================================= */
3855 /* === "Expert" routines =================================================== */
3856 /* ========================================================================= */
3857
3858 /* The following routines are visible outside this routine, but are not meant
3859 * to be called by the user. They are meant for a future version of UMFPACK,
3860 * to replace UMFPACK internal routines with a similar name.
3861 */
3862
3863
3864 /* ========================================================================== */
3865 /* === CCOLAMD_apply_order ================================================== */
3866 /* ========================================================================== */
3867
3868 /*
3869 * Apply post-ordering of supernodal elimination tree.
3870 */
3871
CCOLAMD_apply_order(Int Front[],const Int Order[],Int Temp[],Int nn,Int nfr)3872 GLOBAL void CCOLAMD_apply_order
3873 (
3874 Int Front [ ], /* of size nn on input, size nfr on output */
3875 const Int Order [ ], /* Order [i] = k, i in the range 0..nn-1,
3876 * and k in the range 0..nfr-1, means that node
3877 * i is the kth node in the postordered tree. */
3878 Int Temp [ ], /* workspace of size nfr */
3879 Int nn, /* nodes are numbered in the range 0..nn-1 */
3880 Int nfr /* the number of nodes actually in use */
3881 )
3882 {
3883 Int i, k ;
3884 for (i = 0 ; i < nn ; i++)
3885 {
3886 k = Order [i] ;
3887 ASSERT (k >= EMPTY && k < nfr) ;
3888 if (k != EMPTY)
3889 {
3890 Temp [k] = Front [i] ;
3891 }
3892 }
3893
3894 for (k = 0 ; k < nfr ; k++)
3895 {
3896 Front [k] = Temp [k] ;
3897 }
3898 }
3899
3900
3901 /* ========================================================================== */
3902 /* === CCOLAMD_fsize ======================================================== */
3903 /* ========================================================================== */
3904
3905 /* Determine the largest frontal matrix size for each subtree.
3906 * Only required to sort the children of each
3907 * node prior to postordering the column elimination tree. */
3908
CCOLAMD_fsize(Int nn,Int Fsize[],Int Fnrows[],Int Fncols[],Int Parent[],Int Npiv[])3909 GLOBAL void CCOLAMD_fsize
3910 (
3911 Int nn,
3912 Int Fsize [ ],
3913 Int Fnrows [ ],
3914 Int Fncols [ ],
3915 Int Parent [ ],
3916 Int Npiv [ ]
3917 )
3918 {
3919 double dr, dc ;
3920 Int j, parent, frsize, r, c ;
3921
3922 for (j = 0 ; j < nn ; j++)
3923 {
3924 Fsize [j] = EMPTY ;
3925 }
3926
3927 /* ---------------------------------------------------------------------- */
3928 /* find max front size for tree rooted at node j, for each front j */
3929 /* ---------------------------------------------------------------------- */
3930
3931 DEBUG1 (("\n\n========================================FRONTS:\n")) ;
3932 for (j = 0 ; j < nn ; j++)
3933 {
3934 if (Npiv [j] > 0)
3935 {
3936 /* this is a frontal matrix */
3937 parent = Parent [j] ;
3938 r = Fnrows [j] ;
3939 c = Fncols [j] ;
3940 /* avoid integer overflow */
3941 dr = (double) r ;
3942 dc = (double) c ;
3943 frsize = (INT_OVERFLOW (dr * dc)) ? Int_MAX : (r * c) ;
3944 DEBUG1 ((""ID" : npiv "ID" size "ID" parent "ID" ",
3945 j, Npiv [j], frsize, parent)) ;
3946 Fsize [j] = MAX (Fsize [j], frsize) ;
3947 DEBUG1 (("Fsize [j = "ID"] = "ID"\n", j, Fsize [j])) ;
3948 if (parent != EMPTY)
3949 {
3950 /* find the maximum frontsize of self and children */
3951 ASSERT (Npiv [parent] > 0) ;
3952 ASSERT (parent > j) ;
3953 Fsize [parent] = MAX (Fsize [parent], Fsize [j]) ;
3954 DEBUG1 (("Fsize [parent = "ID"] = "ID"\n",
3955 parent, Fsize [parent]));
3956 }
3957 }
3958 }
3959 DEBUG1 (("fsize done\n")) ;
3960 }
3961
3962
3963 /* ========================================================================= */
3964 /* === CCOLAMD_postorder =================================================== */
3965 /* ========================================================================= */
3966
3967 /* Perform a postordering (via depth-first search) of an assembly tree. */
3968
CCOLAMD_postorder(Int nn,Int Parent[],Int Nv[],Int Fsize[],Int Order[],Int Child[],Int Sibling[],Int Stack[],Int Front_cols[],Int cmember[])3969 GLOBAL void CCOLAMD_postorder
3970 (
3971 /* inputs, not modified on output: */
3972 Int nn, /* nodes are in the range 0..nn-1 */
3973 Int Parent [ ], /* Parent [j] is the parent of j, or EMPTY if root */
3974 Int Nv [ ], /* Nv [j] > 0 number of pivots represented by node j,
3975 * or zero if j is not a node. */
3976 Int Fsize [ ], /* Fsize [j]: size of node j */
3977
3978 /* output, not defined on input: */
3979 Int Order [ ], /* output post-order */
3980
3981 /* workspaces of size nn: */
3982 Int Child [ ],
3983 Int Sibling [ ],
3984 Int Stack [ ],
3985 Int Front_cols [ ],
3986
3987 /* input, not modified on output: */
3988 Int cmember [ ]
3989 )
3990 {
3991 Int i, j, k, parent, frsize, f, fprev, maxfrsize, bigfprev, bigf, fnext ;
3992
3993 for (j = 0 ; j < nn ; j++)
3994 {
3995 Child [j] = EMPTY ;
3996 Sibling [j] = EMPTY ;
3997 }
3998
3999 /* --------------------------------------------------------------------- */
4000 /* place the children in link lists - bigger elements tend to be last */
4001 /* --------------------------------------------------------------------- */
4002
4003 for (j = nn-1 ; j >= 0 ; j--)
4004 {
4005 if (Nv [j] > 0)
4006 {
4007 /* this is an element */
4008 parent = Parent [j] ;
4009 if (parent != EMPTY)
4010 {
4011 /* place the element in link list of the children its parent */
4012 /* bigger elements will tend to be at the end of the list */
4013 Sibling [j] = Child [parent] ;
4014 if (CMEMBER (Front_cols[parent]) == CMEMBER (Front_cols[j]))
4015 {
4016 Child [parent] = j ;
4017 }
4018 }
4019 }
4020 }
4021
4022 #ifndef NDEBUG
4023 {
4024 Int nels, ff, nchild ;
4025 DEBUG1 (("\n\n================================ ccolamd_postorder:\n"));
4026 nels = 0 ;
4027 for (j = 0 ; j < nn ; j++)
4028 {
4029 if (Nv [j] > 0)
4030 {
4031 DEBUG1 ((""ID" : nels "ID" npiv "ID" size "ID
4032 " parent "ID" maxfr "ID"\n", j, nels,
4033 Nv [j], Fsize [j], Parent [j], Fsize [j])) ;
4034 /* this is an element */
4035 /* dump the link list of children */
4036 nchild = 0 ;
4037 DEBUG1 ((" Children: ")) ;
4038 for (ff = Child [j] ; ff != EMPTY ; ff = Sibling [ff])
4039 {
4040 DEBUG1 ((ID" ", ff)) ;
4041 nchild++ ;
4042 ASSERT (nchild < nn) ;
4043 }
4044 DEBUG1 (("\n")) ;
4045 parent = Parent [j] ;
4046 nels++ ;
4047 }
4048 }
4049 }
4050 #endif
4051
4052 /* --------------------------------------------------------------------- */
4053 /* place the largest child last in the list of children for each node */
4054 /* --------------------------------------------------------------------- */
4055
4056 for (i = 0 ; i < nn ; i++)
4057 {
4058 if (Nv [i] > 0 && Child [i] != EMPTY)
4059 {
4060
4061 #ifndef NDEBUG
4062 Int nchild ;
4063 DEBUG1 (("Before partial sort, element "ID"\n", i)) ;
4064 nchild = 0 ;
4065 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4066 {
4067 DEBUG1 ((" f: "ID" size: "ID"\n", f, Fsize [f])) ;
4068 nchild++ ;
4069 }
4070 #endif
4071
4072 /* find the biggest element in the child list */
4073 fprev = EMPTY ;
4074 maxfrsize = EMPTY ;
4075 bigfprev = EMPTY ;
4076 bigf = EMPTY ;
4077 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4078 {
4079 frsize = Fsize [f] ;
4080 if (frsize >= maxfrsize)
4081 {
4082 /* this is the biggest seen so far */
4083 maxfrsize = frsize ;
4084 bigfprev = fprev ;
4085 bigf = f ;
4086 }
4087 fprev = f ;
4088 }
4089
4090 fnext = Sibling [bigf] ;
4091
4092 DEBUG1 (("bigf "ID" maxfrsize "ID" bigfprev "ID" fnext "ID
4093 " fprev " ID"\n", bigf, maxfrsize, bigfprev, fnext, fprev)) ;
4094
4095 if (fnext != EMPTY)
4096 {
4097 /* if fnext is EMPTY then bigf is already at the end of list */
4098
4099 if (bigfprev == EMPTY)
4100 {
4101 /* delete bigf from the element of the list */
4102 Child [i] = fnext ;
4103 }
4104 else
4105 {
4106 /* delete bigf from the middle of the list */
4107 Sibling [bigfprev] = fnext ;
4108 }
4109
4110 /* put bigf at the end of the list */
4111 Sibling [bigf] = EMPTY ;
4112 Sibling [fprev] = bigf ;
4113 }
4114
4115 #ifndef NDEBUG
4116 DEBUG1 (("After partial sort, element "ID"\n", i)) ;
4117 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4118 {
4119 DEBUG1 ((" "ID" "ID"\n", f, Fsize [f])) ;
4120 nchild-- ;
4121 }
4122 #endif
4123 }
4124 }
4125
4126 /* --------------------------------------------------------------------- */
4127 /* postorder the assembly tree */
4128 /* --------------------------------------------------------------------- */
4129
4130 for (i = 0 ; i < nn ; i++)
4131 {
4132 Order [i] = EMPTY ;
4133 }
4134
4135 k = 0 ;
4136
4137 for (i = 0 ; i < nn ; i++)
4138 {
4139 if ((Parent [i] == EMPTY
4140 || (CMEMBER (Front_cols [Parent [i]]) != CMEMBER (Front_cols [i])))
4141 && Nv [i] > 0)
4142 {
4143 DEBUG1 (("Root of assembly tree "ID"\n", i)) ;
4144 k = CCOLAMD_post_tree (i, k, Child, Sibling, Order, Stack) ;
4145 }
4146 }
4147 }
4148
4149
4150 /* ========================================================================= */
4151 /* === CCOLAMD_post_tree =================================================== */
4152 /* ========================================================================= */
4153
4154 /* Post-ordering of a supernodal column elimination tree. */
4155
CCOLAMD_post_tree(Int root,Int k,Int Child[],const Int Sibling[],Int Order[],Int Stack[])4156 GLOBAL Int CCOLAMD_post_tree
4157 (
4158 Int root, /* root of the tree */
4159 Int k, /* start numbering at k */
4160 Int Child [ ], /* input argument of size nn, undefined on
4161 * output. Child [i] is the head of a link
4162 * list of all nodes that are children of node
4163 * i in the tree. */
4164 const Int Sibling [ ], /* input argument of size nn, not modified.
4165 * If f is a node in the link list of the
4166 * children of node i, then Sibling [f] is the
4167 * next child of node i.
4168 */
4169 Int Order [ ], /* output order, of size nn. Order [i] = k
4170 * if node i is the kth node of the reordered
4171 * tree. */
4172 Int Stack [ ] /* workspace of size nn */
4173 )
4174 {
4175 Int f, head, h, i ;
4176
4177 #if 0
4178 /* --------------------------------------------------------------------- */
4179 /* recursive version (Stack [ ] is not used): */
4180 /* --------------------------------------------------------------------- */
4181
4182 /* this is simple, but can cause stack overflow if nn is large */
4183 i = root ;
4184 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4185 {
4186 k = CCOLAMD_post_tree (f, k, Child, Sibling, Order, Stack, nn) ;
4187 }
4188 Order [i] = k++ ;
4189 return (k) ;
4190 #endif
4191
4192 /* --------------------------------------------------------------------- */
4193 /* non-recursive version, using an explicit stack */
4194 /* --------------------------------------------------------------------- */
4195
4196 /* push root on the stack */
4197 head = 0 ;
4198 Stack [0] = root ;
4199
4200 while (head >= 0)
4201 {
4202 /* get head of stack */
4203 i = Stack [head] ;
4204 DEBUG1 (("head of stack "ID" \n", i)) ;
4205
4206 if (Child [i] != EMPTY)
4207 {
4208 /* the children of i are not yet ordered */
4209 /* push each child onto the stack in reverse order */
4210 /* so that small ones at the head of the list get popped first */
4211 /* and the biggest one at the end of the list gets popped last */
4212 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4213 {
4214 head++ ;
4215 }
4216 h = head ;
4217 for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
4218 {
4219 ASSERT (h > 0) ;
4220 Stack [h--] = f ;
4221 DEBUG1 (("push "ID" on stack\n", f)) ;
4222 }
4223 ASSERT (Stack [h] == i) ;
4224
4225 /* delete child list so that i gets ordered next time we see it */
4226 Child [i] = EMPTY ;
4227 }
4228 else
4229 {
4230 /* the children of i (if there were any) are already ordered */
4231 /* remove i from the stack and order it. Front i is kth front */
4232 head-- ;
4233 DEBUG1 (("pop "ID" order "ID"\n", i, k)) ;
4234 Order [i] = k++ ;
4235 }
4236
4237 #ifndef NDEBUG
4238 DEBUG1 (("\nStack:")) ;
4239 for (h = head ; h >= 0 ; h--)
4240 {
4241 Int j = Stack [h] ;
4242 DEBUG1 ((" "ID, j)) ;
4243 }
4244 DEBUG1 (("\n\n")) ;
4245 #endif
4246
4247 }
4248 return (k) ;
4249 }
4250
4251
4252
4253 /* ========================================================================== */
4254 /* === CCOLAMD debugging routines =========================================== */
4255 /* ========================================================================== */
4256
4257 /* When debugging is disabled, the remainder of this file is ignored. */
4258
4259 #ifndef NDEBUG
4260
4261
4262 /* ========================================================================== */
4263 /* === debug_structures ===================================================== */
4264 /* ========================================================================== */
4265
4266 /*
4267 * At this point, all empty rows and columns are dead. All live columns
4268 * are "clean" (containing no dead rows) and simplicial (no supercolumns
4269 * yet). Rows may contain dead columns, but all live rows contain at
4270 * least one live column.
4271 */
4272
debug_structures(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[],Int cmember[],Int cset_start[])4273 PRIVATE void debug_structures
4274 (
4275 /* === Parameters ======================================================= */
4276
4277 Int n_row,
4278 Int n_col,
4279 CColamd_Row Row [ ],
4280 CColamd_Col Col [ ],
4281 Int A [ ],
4282 Int cmember [ ],
4283 Int cset_start [ ]
4284 )
4285 {
4286 /* === Local variables ================================================== */
4287
4288 Int i ;
4289 Int c ;
4290 Int *cp ;
4291 Int *cp_end ;
4292 Int len ;
4293 Int score ;
4294 Int r ;
4295 Int *rp ;
4296 Int *rp_end ;
4297 Int deg ;
4298 Int cs ;
4299
4300 /* === Check A, Row, and Col ============================================ */
4301
4302 for (c = 0 ; c < n_col ; c++)
4303 {
4304 if (COL_IS_ALIVE (c))
4305 {
4306 len = Col [c].length ;
4307 score = Col [c].shared2.score ;
4308 DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
4309 ASSERT (len > 0) ;
4310 ASSERT (score >= 0) ;
4311 ASSERT (Col [c].shared1.thickness == 1) ;
4312 cp = &A [Col [c].start] ;
4313 cp_end = cp + len ;
4314 while (cp < cp_end)
4315 {
4316 r = *cp++ ;
4317 ASSERT (ROW_IS_ALIVE (r)) ;
4318 }
4319 }
4320 else
4321 {
4322 i = Col [c].shared2.order ;
4323 cs = CMEMBER (c) ;
4324 ASSERT (i >= cset_start [cs] && i < cset_start [cs+1]) ;
4325 }
4326 }
4327
4328 for (r = 0 ; r < n_row ; r++)
4329 {
4330 if (ROW_IS_ALIVE (r))
4331 {
4332 i = 0 ;
4333 len = Row [r].length ;
4334 deg = Row [r].shared1.degree ;
4335 ASSERT (len > 0) ;
4336 ASSERT (deg > 0) ;
4337 rp = &A [Row [r].start] ;
4338 rp_end = rp + len ;
4339 while (rp < rp_end)
4340 {
4341 c = *rp++ ;
4342 if (COL_IS_ALIVE (c))
4343 {
4344 i++ ;
4345 }
4346 }
4347 ASSERT (i > 0) ;
4348 }
4349 }
4350 }
4351
4352
4353 /* ========================================================================== */
4354 /* === debug_deg_lists ====================================================== */
4355 /* ========================================================================== */
4356
4357 /*
4358 * Prints the contents of the degree lists. Counts the number of columns
4359 * in the degree list and compares it to the total it should have. Also
4360 * checks the row degrees.
4361 */
4362
debug_deg_lists(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int head[],Int min_score,Int should,Int max_deg)4363 PRIVATE void debug_deg_lists
4364 (
4365 /* === Parameters ======================================================= */
4366
4367 Int n_row,
4368 Int n_col,
4369 CColamd_Row Row [ ],
4370 CColamd_Col Col [ ],
4371 Int head [ ],
4372 Int min_score,
4373 Int should,
4374 Int max_deg
4375 )
4376
4377 {
4378 /* === Local variables ================================================== */
4379
4380 Int deg ;
4381 Int col ;
4382 Int have ;
4383 Int row ;
4384
4385 /* === Check the degree lists =========================================== */
4386
4387 if (n_col > 10000 && ccolamd_debug <= 0)
4388 {
4389 return ;
4390 }
4391 have = 0 ;
4392 DEBUG4 (("Degree lists: "ID"\n", min_score)) ;
4393 for (deg = 0 ; deg <= n_col ; deg++)
4394 {
4395 col = head [deg] ;
4396 if (col == EMPTY)
4397 {
4398 continue ;
4399 }
4400 DEBUG4 (("%d:", deg)) ;
4401 ASSERT (Col [col].shared3.prev == EMPTY) ;
4402 while (col != EMPTY)
4403 {
4404 DEBUG4 ((" "ID"", col)) ;
4405 have += Col [col].shared1.thickness ;
4406 ASSERT (COL_IS_ALIVE (col)) ;
4407 col = Col [col].shared4.degree_next ;
4408 }
4409 DEBUG4 (("\n")) ;
4410 }
4411 DEBUG4 (("should "ID" have "ID"\n", should, have)) ;
4412 ASSERT (should == have) ;
4413
4414 /* === Check the row degrees ============================================ */
4415
4416 if (n_row > 10000 && ccolamd_debug <= 0)
4417 {
4418 return ;
4419 }
4420 for (row = 0 ; row < n_row ; row++)
4421 {
4422 if (ROW_IS_ALIVE (row))
4423 {
4424 ASSERT (Row [row].shared1.degree <= max_deg) ;
4425 }
4426 }
4427 }
4428
4429
4430 /* ========================================================================== */
4431 /* === debug_mark =========================================================== */
4432 /* ========================================================================== */
4433
4434 /*
4435 * Ensures that the tag_mark is less that the maximum and also ensures that
4436 * each entry in the mark array is less than the tag mark.
4437 */
4438
debug_mark(Int n_row,CColamd_Row Row[],Int tag_mark,Int max_mark)4439 PRIVATE void debug_mark
4440 (
4441 /* === Parameters ======================================================= */
4442
4443 Int n_row,
4444 CColamd_Row Row [ ],
4445 Int tag_mark,
4446 Int max_mark
4447 )
4448 {
4449 /* === Local variables ================================================== */
4450
4451 Int r ;
4452
4453 /* === Check the Row marks ============================================== */
4454
4455 ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
4456 if (n_row > 10000 && ccolamd_debug <= 0)
4457 {
4458 return ;
4459 }
4460 for (r = 0 ; r < n_row ; r++)
4461 {
4462 ASSERT (Row [r].shared2.mark < tag_mark) ;
4463 }
4464 }
4465
4466
4467 /* ========================================================================== */
4468 /* === debug_matrix ========================================================= */
4469 /* ========================================================================== */
4470
4471 /* Prints out the contents of the columns and the rows. */
4472
debug_matrix(Int n_row,Int n_col,CColamd_Row Row[],CColamd_Col Col[],Int A[])4473 PRIVATE void debug_matrix
4474 (
4475 /* === Parameters ======================================================= */
4476
4477 Int n_row,
4478 Int n_col,
4479 CColamd_Row Row [ ],
4480 CColamd_Col Col [ ],
4481 Int A [ ]
4482 )
4483 {
4484 /* === Local variables ================================================== */
4485
4486 Int r ;
4487 Int c ;
4488 Int *rp ;
4489 Int *rp_end ;
4490 Int *cp ;
4491 Int *cp_end ;
4492
4493 /* === Dump the rows and columns of the matrix ========================== */
4494
4495 if (ccolamd_debug < 3)
4496 {
4497 return ;
4498 }
4499 DEBUG3 (("DUMP MATRIX:\n")) ;
4500 for (r = 0 ; r < n_row ; r++)
4501 {
4502 DEBUG3 (("Row "ID" alive? "ID"\n", r, ROW_IS_ALIVE (r))) ;
4503 if (ROW_IS_DEAD (r))
4504 {
4505 continue ;
4506 }
4507
4508 DEBUG3 (("start "ID" length "ID" degree "ID"\nthickness "ID"\n",
4509 Row [r].start, Row [r].length, Row [r].shared1.degree,
4510 Row [r].thickness)) ;
4511
4512 rp = &A [Row [r].start] ;
4513 rp_end = rp + Row [r].length ;
4514 while (rp < rp_end)
4515 {
4516 c = *rp++ ;
4517 DEBUG4 ((" "ID" col "ID"\n", COL_IS_ALIVE (c), c)) ;
4518 }
4519 }
4520
4521 for (c = 0 ; c < n_col ; c++)
4522 {
4523 DEBUG3 (("Col "ID" alive? "ID"\n", c, COL_IS_ALIVE (c))) ;
4524 if (COL_IS_DEAD (c))
4525 {
4526 continue ;
4527 }
4528 DEBUG3 (("start "ID" length "ID" shared1 "ID" shared2 "ID"\n",
4529 Col [c].start, Col [c].length,
4530 Col [c].shared1.thickness, Col [c].shared2.score)) ;
4531 cp = &A [Col [c].start] ;
4532 cp_end = cp + Col [c].length ;
4533 while (cp < cp_end)
4534 {
4535 r = *cp++ ;
4536 DEBUG4 ((" "ID" row "ID"\n", ROW_IS_ALIVE (r), r)) ;
4537 }
4538 }
4539 }
4540
4541
4542 /* ========================================================================== */
4543 /* === dump_super =========================================================== */
4544 /* ========================================================================== */
4545
dump_super(Int super_c,CColamd_Col Col[],Int n_col)4546 PRIVATE void dump_super
4547 (
4548 Int super_c,
4549 CColamd_Col Col [ ],
4550 Int n_col
4551 )
4552 {
4553 Int col, ncols ;
4554
4555 DEBUG1 ((" =[ ")) ;
4556 ncols = 0 ;
4557 for (col = super_c ; col != EMPTY ; col = Col [col].nextcol)
4558 {
4559 DEBUG1 ((" "ID, col)) ;
4560 ASSERT (col >= 0 && col < n_col) ;
4561 if (col != super_c)
4562 {
4563 ASSERT (COL_IS_DEAD (col)) ;
4564 }
4565 if (Col [col].nextcol == EMPTY)
4566 {
4567 ASSERT (col == Col [super_c].lastcol) ;
4568 }
4569 ncols++ ;
4570 ASSERT (ncols <= Col [super_c].shared1.thickness) ;
4571 }
4572 ASSERT (ncols == Col [super_c].shared1.thickness) ;
4573 DEBUG1 (("]\n")) ;
4574 }
4575
4576
4577 /* ========================================================================== */
4578 /* === ccolamd_get_debug ==================================================== */
4579 /* ========================================================================== */
4580
ccolamd_get_debug(char * method)4581 PRIVATE void ccolamd_get_debug
4582 (
4583 char *method
4584 )
4585 {
4586 FILE *debug_file ;
4587 ccolamd_debug = 0 ; /* no debug printing */
4588
4589 /* Read debug info from the debug file. */
4590 debug_file = fopen ("debug", "r") ;
4591 if (debug_file)
4592 {
4593 (void) fscanf (debug_file, ""ID"", &ccolamd_debug) ;
4594 (void) fclose (debug_file) ;
4595 }
4596
4597 DEBUG0 ((":")) ;
4598 DEBUG1 (("%s: debug version, D = "ID" (THIS WILL BE SLOW!)\n",
4599 method, ccolamd_debug)) ;
4600 DEBUG1 ((" Debug printing level: "ID"\n", ccolamd_debug)) ;
4601 }
4602
4603 #endif
4604