ccolamd.c
Go to the documentation of this file.
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  {
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 
881 (
882  Int n_row,
883  Int n_col,
884  CColamd_Row Row [ ],
885  CColamd_Col Col [ ],
886  Int A [ ],
887  Int p [ ],
889 ) ;
890 
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 
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 
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 
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 
971 (
972  Int tag_mark,
973  Int max_mark,
974  Int n_row,
975  CColamd_Row Row [ ]
976 ) ;
977 
979 (
980  char *method,
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 */
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 */
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 
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 
1083 PUBLIC size_t CCOLAMD_recommended /* returns recommended value of Alen. */
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 
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 
1144 PUBLIC Int CSYMAMD_MAIN /* return TRUE if OK, FALSE otherwise */
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  {
1210  DEBUG1 (("csymamd: A not present\n")) ;
1211  return (FALSE) ;
1212  }
1213 
1214  if (!p) /* p is not present */
1215  {
1217  DEBUG1 (("csymamd: p not present\n")) ;
1218  return (FALSE) ;
1219  }
1220 
1221  if (n < 0) /* n must be >= 0 */
1222  {
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  {
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  {
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  {
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  {
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 */
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 */
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. */
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  {
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 
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 
1505 PUBLIC Int CCOLAMD_2 /* returns TRUE if successful, FALSE otherwise */
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  {
1589  DEBUG1 (("ccolamd: A not present\n")) ;
1590  return (FALSE) ;
1591  }
1592 
1593  if (!p) /* p is not present */
1594  {
1596  DEBUG1 (("ccolamd: p not present\n")) ;
1597  return (FALSE) ;
1598  }
1599 
1600  if (n_row < 0) /* n_row must be >= 0 */
1601  {
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  {
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  {
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  {
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 */
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  {
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 
1974  Int stats [CCOLAMD_STATS]
1975 )
1976 {
1977  print_report ("ccolamd", stats) ;
1978 }
1979 
1980 
1981 /* ========================================================================== */
1982 /* === symamd_report ======================================================== */
1983 /* ========================================================================== */
1984 
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 
2014 PRIVATE Int init_rows_cols /* returns TRUE if OK, or FALSE otherwise */
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 */
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  {
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. */
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 
2259 PRIVATE void init_scoring
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 
2545 PRIVATE Int find_ordering /* return the number of garbage collections */
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 
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 
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 ;
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 
3545 PRIVATE Int garbage_collection /* returns the new value of pfree */
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 
3677 PRIVATE Int clear_mark /* return the new value for tag_mark */
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 
3714 PRIVATE void print_report
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,
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 
3749 
3751  "Matrix has unsorted or duplicate row indices.\n")) ;
3752 
3754  "%s: duplicate or out-of-order row indices: "ID"\n",
3755  method, i3)) ;
3756 
3758  "%s: last seen duplicate or out-of-order row: "ID"\n",
3759  method, INDEX (i2))) ;
3760 
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 
3772  "%s: number of dense or empty rows ignored: "ID"\n",
3773  method, stats [CCOLAMD_DENSE_ROW])) ;
3774 
3776  "%s: number of dense or empty columns ignored: "ID"\n",
3777  method, stats [CCOLAMD_DENSE_COL])) ;
3778 
3780  "%s: number of garbage collections performed: "ID"\n",
3781  method, stats [CCOLAMD_DEFRAG_COUNT])) ;
3782  break ;
3783 
3785 
3787  "Array A (row indices of matrix) not present.\n")) ;
3788  break ;
3789 
3791 
3793  "Array p (column pointers for matrix) not present.\n")) ;
3794  break ;
3795 
3797 
3798  SUITESPARSE_PRINTF(("Invalid number of rows ("ID").\n", i1)) ;
3799  break ;
3800 
3802 
3803  SUITESPARSE_PRINTF(("Invalid number of columns ("ID").\n", i1)) ;
3804  break ;
3805 
3807 
3809  "Invalid number of nonzero entries ("ID").\n", i1)) ;
3810  break ;
3811 
3813 
3815  "Invalid column pointer, p [0] = "ID", must be 0.\n", i1)) ;
3816  break ;
3817 
3819 
3820  SUITESPARSE_PRINTF(("Array A too small.\n")) ;
3822  " Need Alen >= "ID", but given only Alen = "ID".\n",
3823  i1, i2)) ;
3824  break ;
3825 
3827 
3829  "Column "ID" has a negative number of entries ("ID").\n",
3830  INDEX (i1), i2)) ;
3831  break ;
3832 
3834 
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 
3842 
3843  SUITESPARSE_PRINTF(("Out of memory.\n")) ;
3844  break ;
3845 
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 
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 
3909 GLOBAL void CCOLAMD_fsize
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 
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 
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 
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 
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 
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 
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 
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 
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
#define CCOLAMD_report
Definition: ccolamd.c:649
#define CCOLAMD_ERROR_row_index_out_of_bounds
Definition: ccolamd.h:105
#define CCOLAMD_SUB_VERSION
Definition: ccolamd.h:47
union CColamd_Col_struct::@0 shared1
PRIVATE void 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)
Definition: ccolamd.c:2260
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51
#define CMEMBER(c)
Definition: ccolamd.c:738
static char lower
#define CCOLAMD_ERROR_p0_nonzero
Definition: ccolamd.h:102
#define CCOLAMD_ERROR_col_length_negative
Definition: ccolamd.h:104
Scalar * b
Definition: benchVecAdd.cpp:17
#define CCOLAMD_INFO2
Definition: ccolamd.h:82
struct CColamd_Row_struct CColamd_Row
#define KILL_NON_PRINCIPAL_COL(c)
Definition: ccolamd.c:770
#define CCOLAMD_post_tree
Definition: ccolamd.c:646
#define CCOLAMD_ERROR_p_not_present
Definition: ccolamd.h:98
#define INT_OVERFLOW(x)
Definition: ccolamd.c:744
#define ROW_IS_ALIVE(r)
Definition: ccolamd.c:764
#define CCOLAMD_DENSE_ROW
Definition: ccolamd.h:63
#define Int_MAX
Definition: ccolamd.c:638
#define FALSE
Definition: ccolamd.c:751
PRIVATE Int init_rows_cols(Int n_row, Int n_col, CColamd_Row Row [], CColamd_Col Col [], Int A [], Int p [], Int stats [CCOLAMD_STATS])
Definition: ccolamd.c:2015
int n
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
#define CCOLAMD_ERROR_nnz_negative
Definition: ccolamd.h:101
#define CCOLAMD_ERROR_A_not_present
Definition: ccolamd.h:97
PRIVATE void detect_super_cols(CColamd_Col Col [], Int A [], Int head [], Int row_start, Int row_length, Int in_set [])
Definition: ccolamd.c:3374
#define DEBUG4(params)
Definition: ccolamd.c:870
#define DENSE_DEGREE(alpha, n)
Definition: ccolamd.c:735
#define CCOLAMD_recommended
Definition: ccolamd.c:640
#define CSYMAMD_report
Definition: ccolamd.c:650
#define CCOLAMD_MAIN_VERSION
Definition: ccolamd.h:46
#define MIN(a, b)
Definition: ccolamd.c:728
#define GLOBAL
Definition: ccolamd.c:731
#define CSYMAMD_MAIN
Definition: ccolamd.c:648
if n return
#define CCOLAMD_NEWLY_EMPTY_COL
Definition: ccolamd.h:92
#define EMPTY
Definition: ccolamd.c:726
#define CCOLAMD_ERROR_nrow_negative
Definition: ccolamd.h:99
#define CCOLAMD_OK
Definition: ccolamd.h:95
#define ROW_IS_MARKED_DEAD(row_mark)
Definition: ccolamd.c:763
#define CCOLAMD_ERROR_invalid_cmember
Definition: ccolamd.h:107
#define CCOLAMD_ERROR_ncol_negative
Definition: ccolamd.h:100
#define CCOLAMD_fsize
Definition: ccolamd.c:647
bool stats
#define CCOLAMD_LU
Definition: ccolamd.h:72
#define DEBUG1(params)
Definition: ccolamd.c:867
#define CCOLAMD_EMPTY_COL
Definition: ccolamd.h:88
struct CColamd_Col_struct CColamd_Col
#define ASSERT(expression)
Definition: ccolamd.c:872
#define CCOLAMD_NEWLY_EMPTY_ROW
Definition: ccolamd.h:90
union CColamd_Col_struct::@3 shared4
union CColamd_Col_struct::@1 shared2
#define DEBUG0(params)
Definition: ccolamd.c:866
#define CCOLAMD_apply_order
Definition: ccolamd.c:644
#define PRIVATE
Definition: ccolamd.c:733
#define CCOLAMD_set_defaults
Definition: ccolamd.c:641
m row(1)
#define CCOLAMD_ERROR_out_of_memory
Definition: ccolamd.h:106
#define CCOLAMD_STATUS
Definition: ccolamd.h:78
#define CCOLAMD_DATE
Definition: ccolamd.h:44
static size_t ccolamd_need(Int nnz, Int n_row, Int n_col, int *ok)
Definition: ccolamd.c:1047
static size_t t_add(size_t a, size_t b, int *ok)
Definition: ccolamd.c:1016
PRIVATE Int clear_mark(Int tag_mark, Int max_mark, Int n_row, CColamd_Row Row [])
Definition: ccolamd.c:3678
PRIVATE Int find_ordering(Int n_row, Int n_col, Int Alen, CColamd_Row Row [], CColamd_Col Col [], Int A [], Int head [], Int max_deg, Int pfree, Int cset [], Int cset_start [], 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)
Definition: ccolamd.c:2546
#define ONES_COMPLEMENT(r)
Definition: ccolamd.c:747
idx_t idx_t idx_t idx_t idx_t * perm
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
RealScalar s
#define CCOLAMD_MAIN
Definition: ccolamd.c:643
PRIVATE void print_report(char *method, Int stats [CCOLAMD_STATS])
Definition: ccolamd.c:3715
#define MAX(a, b)
Definition: ccolamd.c:727
#define NULL
Definition: ccolamd.c:609
std::vector< double > Row
#define NDEBUG
Definition: ccolamd.c:582
#define CCOLAMD_DENSE_COL
Definition: ccolamd.h:66
#define TRUE
Definition: ccolamd.c:750
#define DEBUG3(params)
Definition: ccolamd.c:869
union CColamd_Col_struct::@2 shared3
const double h
#define CCOLAMD_EMPTY_ROW
Definition: ccolamd.h:86
idx_t * nn
#define INDEX(i)
Definition: ccolamd.c:782
#define DEBUG2(params)
Definition: ccolamd.c:868
#define CCOLAMD_AGGRESSIVE
Definition: ccolamd.h:69
static size_t t_mult(size_t a, size_t k, int *ok)
Definition: ccolamd.c:1023
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE FixedSegmentReturnType< internal::get_fixed_value< NType >::value >::Type head(NType n)
#define CCOLAMD_ERROR_A_too_small
Definition: ccolamd.h:103
float * p
#define CCOLAMD_KNOBS
Definition: ccolamd.h:57
#define KILL_ROW(r)
Definition: ccolamd.c:768
#define COL_IS_ALIVE(c)
Definition: ccolamd.c:766
#define ID
Definition: ccolamd.c:637
#define CCOLAMD_postorder
Definition: ccolamd.c:645
#define CCOLAMD_C(n_col, ok)
Definition: ccolamd.c:1034
m col(1)
union CColamd_Row_struct::@4 shared1
#define CCOLAMD_OK_BUT_JUMBLED
Definition: ccolamd.h:96
PRIVATE Int garbage_collection(Int n_row, Int n_col, CColamd_Row Row [], CColamd_Col Col [], Int A [], Int *pfree)
Definition: ccolamd.c:3546
#define ROW_IS_DEAD(r)
Definition: ccolamd.c:762
#define CCOLAMD_R(n_row, ok)
Definition: ccolamd.c:1037
#define CCOLAMD_INFO3
Definition: ccolamd.h:83
size_t len(handle h)
Get the length of a Python object.
Definition: pytypes.h:2244
#define CCOLAMD_2
Definition: ccolamd.c:642
#define SUITESPARSE_PRINTF(params)
#define COL_IS_DEAD(c)
Definition: ccolamd.c:765
#define CCOLAMD_STATS
Definition: ccolamd.h:60
#define KILL_PRINCIPAL_COL(c)
Definition: ccolamd.c:769
std::ptrdiff_t j
#define CCOLAMD_DEFRAG_COUNT
Definition: ccolamd.h:75
Point2 t(10, 10)
#define Int
Definition: ccolamd.c:636
#define PUBLIC
Definition: ccolamd.c:732
#define CCOLAMD_INFO1
Definition: ccolamd.h:81


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:01