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