source:glpk-cmake/src/colamd/colamd.c@1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

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