Amesos Package Browser (Single Doxygen Collection) Development
amesos_ccolamd.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === CCOLAMD/CSYMAMD - a constrained column ordering algorithm ============ */
00003 /* ========================================================================== */
00004
00005 /* ----------------------------------------------------------------------------
00006  * CCOLAMD, Copyright (C) Univ. of Florida.  Authors: Timothy A. Davis,
00007  * Sivasankaran Rajamanickam, and Stefan Larimore
00008  * See License.txt for the Version 2.1 of the GNU Lesser General Public License
00009  * http://www.cise.ufl.edu/research/sparse
00010  * -------------------------------------------------------------------------- */
00011
00012 /*
00013  *  ccolamd:  a constrained approximate minimum degree column ordering
00014  *  algorithm, LU factorization of symmetric or unsymmetric matrices,
00015  *  QR factorization, least squares, interior point methods for
00016  *  linear programming problems, and other related problems.
00017  *
00018  *  csymamd:  a constrained approximate minimum degree ordering algorithm for
00019  *  Cholesky factorization of symmetric matrices.
00020  *
00021  *  Purpose:
00022  *
00023  *  CCOLAMD computes a permutation Q such that the Cholesky factorization of
00024  *  (AQ)'(AQ) has less fill-in and requires fewer floating point operations
00025  *  than A'A.  This also provides a good ordering for sparse partial
00026  *  pivoting methods, P(AQ) = LU, where Q is computed prior to numerical
00027  *  factorization, and P is computed during numerical factorization via
00028  *  conventional partial pivoting with row interchanges.  CCOLAMD is an
00029  *  extension of COLAMD, available as built-in function in MATLAB Version 6,
00030  *  available from MathWorks, Inc. (http://www.mathworks.com).  This
00031  *  routine can be used in place of COLAMD in MATLAB.
00032  *
00033  *  CSYMAMD computes a permutation P of a symmetric matrix A such that the
00034  *  Cholesky factorization of PAP' has less fill-in and requires fewer
00035  *  floating point operations than A.  CSYMAMD constructs a matrix M such
00036  *  that M'M has the same nonzero pattern of A, and then orders the columns
00037  *  of M using colmmd.  The column ordering of M is then returned as the
00038  *  row and column ordering P of A.  CSYMAMD is an extension of SYMAMD.
00039  *
00040  *  Authors:
00041  *
00042  *  Timothy A. Davis and S. Rajamanickam wrote CCOLAMD, based directly on
00043  *  COLAMD by Stefan I. Larimore and Timothy A. Davis, University of
00044  *  Florida.  The algorithm was developed in collaboration with John
00045  *  Gilbert, (UCSB, then at Xerox PARC), and Esmond Ng, (Lawrence Berkeley
00046  *  National Lab, then at Oak Ridge National Laboratory).
00047  *
00048  *  Acknowledgements:
00049  *
00050  *  This work was supported by the National Science Foundation, under
00051  *  grants DMS-9504974 and DMS-9803599, CCR-0203270, and a grant from the
00052  *  Sandia National Laboratory (Dept. of Energy).
00053  *
00054  *  Copyright and License:
00055  *
00056  *  Copyright (c) 1998-2005 by the University of Florida.
00057  *  All Rights Reserved.
00058  *  COLAMD is also available under alternate licenses, contact T. Davis
00059  *  for details.
00060  *
00061  *  This library is free software; you can redistribute it and/or
00062  *  modify it under the terms of the GNU Lesser General Public
00063  *  License as published by the Free Software Foundation; either
00064  *  version 2.1 of the License, or (at your option) any later version.
00065  *
00066  *  This library is distributed in the hope that it will be useful,
00067  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00068  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00069  *  Lesser General Public License for more details.
00070  *
00071  *  You should have received a copy of the GNU Lesser General Public
00072  *  License along with this library; if not, write to the Free Software
00073  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
00074  *  USA
00075  *
00076  *  Permission is hereby granted to use or copy this program under the
00077  *  terms of the GNU LGPL, provided that the Copyright, this License,
00078  *  and the Availability of the original version is retained on all copies.
00079  *  User documentation of any code that uses this code or any modified
00080  *  version of this code must cite the Copyright, this License, the
00081  *  Availability note, and "Used by permission." Permission to modify
00082  *  the code and to distribute modified code is granted, provided the
00083  *  Copyright, this License, and the Availability note are retained,
00084  *  and a notice that the code was modified is included.
00085  *
00086  *  Availability:
00087  *
00088  *  The CCOLAMD/CSYMAMD library is available at
00089  *
00090  *      http://www.cise.ufl.edu/research/sparse/ccolamd/
00091  *
00092  *  This is the http://www.cise.ufl.edu/research/sparse/ccolamd/ccolamd.c
00093  *  file.
00094  *
00095  *   See the ChangeLog file for changes since Version 1.0.
00096  */
00097
00098 /* ========================================================================== */
00099 /* === Description of user-callable routines ================================ */
00100 /* ========================================================================== */
00101
00102 /* CCOLAMD includes both int and UF_long versions of all its routines.  The
00103  * description below is for the int version.   For UF_long, all int arguments
00104  * become UF_long integers.  UF_long is normally defined as long, except for
00105  * WIN64 */
00106
00107 /*  ----------------------------------------------------------------------------
00108  *  ccolamd_recommended:
00109  *  ----------------------------------------------------------------------------
00110  *
00111  *  C syntax:
00112  *
00113  *      #include "ccolamd.h"
00114  *      size_t ccolamd_recommended (int nnz, int n_row, int n_col) ;
00115  *      size_t ccolamd_l_recommended (UF_long nnz, UF_long n_row,
00116  *    UF_long n_col) ;
00117  *
00118  *  Purpose:
00119  *
00120  *      Returns recommended value of Alen for use by ccolamd.  Returns 0
00121  *      if any input argument is negative.  The use of this routine
00122  *      is optional.  Not needed for csymamd, which dynamically allocates
00123  *      its own memory.
00124  *
00125  *  Arguments (all input arguments):
00126  *
00127  *      int nnz ;   Number of nonzeros in the matrix A.  This must
00128  *        be the same value as p [n_col] in the call to
00129  *        ccolamd - otherwise you will get a wrong value
00130  *        of the recommended memory to use.
00131  *
00132  *      int n_row ;   Number of rows in the matrix A.
00133  *
00134  *      int n_col ;   Number of columns in the matrix A.
00135  *
00136  *  ----------------------------------------------------------------------------
00137  *  ccolamd_set_defaults:
00138  *  ----------------------------------------------------------------------------
00139  *
00140  *  C syntax:
00141  *
00142  *      #include "ccolamd.h"
00143  *      ccolamd_set_defaults (double knobs [CCOLAMD_KNOBS]) ;
00144  *      ccolamd_l_set_defaults (double knobs [CCOLAMD_KNOBS]) ;
00145  *
00146  *  Purpose:
00147  *
00148  *      Sets the default parameters.  The use of this routine is optional.
00149  *      Passing a (double *) NULL pointer for the knobs results in the
00150  *      default parameter settings.
00151  *
00152  *  Arguments:
00153  *
00154  *      double knobs [CCOLAMD_KNOBS] ;  Output only.
00155  *
00156  *      knobs [0] and knobs [1] behave differently than they did in COLAMD.
00157  *      The other knobs are new to CCOLAMD.
00158  *
00159  *      knobs [0]: dense row control
00160  *
00161  *    For CCOLAMD, rows with more than
00162  *    max (16, knobs [CCOLAMD_DENSE_ROW] * sqrt (n_col))
00163  *    entries are removed prior to ordering.
00164  *
00165  *    For CSYMAMD, rows and columns with more than
00166  *    max (16, knobs [CCOLAMD_DENSE_ROW] * sqrt (n))
00167  *    entries are removed prior to ordering, and placed last in the
00168  *    output ordering (subject to the constraints).
00169  *
00170  *    If negative, only completely dense rows are removed.  If you
00171  *    intend to use CCOLAMD for a Cholesky factorization of A*A', set
00172  *    knobs [CCOLAMD_DENSE_ROW] to -1, which is more appropriate for
00173  *    that case.
00174  *
00175  *    Default: 10.
00176  *
00177  *      knobs [1]: dense column control
00178  *
00179  *    For CCOLAMD, columns with more than
00180  *    max (16, knobs [CCOLAMD_DENSE_COL] * sqrt (MIN (n_row,n_col)))
00181  *    entries are removed prior to ordering, and placed last in the
00182  *    output column ordering (subject to the constraints).
00183  *    Not used by CSYMAMD.  If negative, only completely dense
00184  *    columns are removed.  Default: 10.
00185  *
00186  *      knobs [2]: aggressive absorption
00187  *
00188  *          knobs [CCOLAMD_AGGRESSIVE] controls whether or not to do
00189  *          aggressive absorption during the ordering.  Default is TRUE
00190  *          (nonzero).  If zero, no aggressive absorption is performed.
00191  *
00192  *      knobs [3]: optimize ordering for LU or Cholesky
00193  *
00194  *    knobs [CCOLAMD_LU] controls an option that optimizes the
00195  *    ordering for the LU of A or the Cholesky factorization of A'A.
00196  *    If TRUE (nonzero), an ordering optimized for LU is performed.
00197  *    If FALSE (zero), an ordering for Cholesky is performed.
00198  *    Default is FALSE.  CSYMAMD ignores this parameter; it always
00199  *    orders for Cholesky.
00200  *
00201  *  ----------------------------------------------------------------------------
00202  *  ccolamd:
00203  *  ----------------------------------------------------------------------------
00204  *
00205  *  C syntax:
00206  *
00207  *      #include "ccolamd.h"
00208  *      int ccolamd (int n_row, int n_col, int Alen, int *A, int *p,
00209  *        double knobs [CCOLAMD_KNOBS], int stats [CCOLAMD_STATS],
00210  *    int *cmember) ;
00211  *
00212  *      UF_long ccolamd_l (UF_long n_row, UF_long n_col, UF_long Alen,
00213  *    UF_long *A, UF_long *p, double knobs [CCOLAMD_KNOBS],
00214  *    UF_long stats [CCOLAMD_STATS], UF_long *cmember) ;
00215  *
00216  *  Purpose:
00217  *
00218  *      Computes a column ordering (Q) of A such that P(AQ)=LU or
00219  *      (AQ)'AQ=LL' have less fill-in and require fewer floating point
00220  *      operations than factorizing the unpermuted matrix A or A'A,
00221  *      respectively.
00222  *
00223  *  Returns:
00224  *
00225  *      TRUE (1) if successful, FALSE (0) otherwise.
00226  *
00227  *  Arguments (for int version):
00228  *
00229  *      int n_row ;   Input argument.
00230  *
00231  *    Number of rows in the matrix A.
00232  *    Restriction:  n_row >= 0.
00233  *    ccolamd returns FALSE if n_row is negative.
00234  *
00235  *      int n_col ;   Input argument.
00236  *
00237  *    Number of columns in the matrix A.
00238  *    Restriction:  n_col >= 0.
00239  *    ccolamd returns FALSE if n_col is negative.
00240  *
00241  *      int Alen ;    Input argument.
00242  *
00243  *    Restriction (see note):
00244  *    Alen >= MAX (2*nnz, 4*n_col) + 17*n_col + 7*n_row + 7, where
00245  *    nnz = p [n_col].  ccolamd returns FALSE if this condition is
00246  *    not met. We recommend about nnz/5 more space for better
00247  *    efficiency.  This restriction makes an modest assumption
00248  *    regarding the size of two typedef'd structures in ccolamd.h.
00249  *    We do, however, guarantee that
00250  *
00251  *        Alen >= ccolamd_recommended (nnz, n_row, n_col)
00252  *
00253  *    will work efficiently.
00254  *
00255  *      int A [Alen] ;  Input argument, undefined on output.
00256  *
00257  *    A is an integer array of size Alen.  Alen must be at least as
00258  *    large as the bare minimum value given above, but this is very
00259  *    low, and can result in excessive run time.  For best
00260  *    performance, we recommend that Alen be greater than or equal to
00261  *    ccolamd_recommended (nnz, n_row, n_col), which adds
00262  *    nnz/5 to the bare minimum value given above.
00263  *
00264  *    On input, the row indices of the entries in column c of the
00265  *    matrix are held in A [(p [c]) ... (p [c+1]-1)].  The row indices
00266  *    in a given column c need not be in ascending order, and
00267  *    duplicate row indices may be be present.  However, ccolamd will
00268  *    work a little faster if both of these conditions are met
00269  *    (ccolamd puts the matrix into this format, if it finds that the
00270  *    the conditions are not met).
00271  *
00272  *    The matrix is 0-based.  That is, rows are in the range 0 to
00273  *    n_row-1, and columns are in the range 0 to n_col-1.  ccolamd
00274  *    returns FALSE if any row index is out of range.
00275  *
00276  *    The contents of A are modified during ordering, and are
00277  *    undefined on output.
00278  *
00279  *      int p [n_col+1] ; Both input and output argument.
00280  *
00281  *    p is an integer array of size n_col+1.  On input, it holds the
00282  *    "pointers" for the column form of the matrix A.  Column c of
00283  *    the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
00284  *    entry, p [0], must be zero, and p [c] <= p [c+1] must hold
00285  *    for all c in the range 0 to n_col-1.  The value nnz = p [n_col]
00286  *    is thus the total number of entries in the pattern of the
00287  *    matrix A.  ccolamd returns FALSE if these conditions are not
00288  *    met.
00289  *
00290  *    On output, if ccolamd returns TRUE, the array p holds the column
00291  *    permutation (Q, for P(AQ)=LU or (AQ)'(AQ)=LL'), where p [0] is
00292  *    the first column index in the new ordering, and p [n_col-1] is
00293  *    the last.  That is, p [k] = j means that column j of A is the
00294  *    kth pivot column, in AQ, where k is in the range 0 to n_col-1
00295  *    (p [0] = j means that column j of A is the first column in AQ).
00296  *
00297  *    If ccolamd returns FALSE, then no permutation is returned, and
00298  *    p is undefined on output.
00299  *
00300  *      double knobs [CCOLAMD_KNOBS] ;  Input argument.
00301  *
00302  *    See ccolamd_set_defaults for a description.
00303  *
00304  *      int stats [CCOLAMD_STATS] ;   Output argument.
00305  *
00306  *    Statistics on the ordering, and error status.
00307  *    See ccolamd.h for related definitions.
00308  *    ccolamd returns FALSE if stats is not present.
00309  *
00310  *    stats [0]:  number of dense or empty rows ignored.
00311  *
00312  *    stats [1]:  number of dense or empty columns ignored (and
00313  *        ordered last in the output permutation p, subject to the
00314  *        constraints).  Note that a row can become "empty" if it
00315  *        contains only "dense" and/or "empty" columns, and similarly
00316  *        a column can become "empty" if it only contains "dense"
00317  *        and/or "empty" rows.
00318  *
00319  *    stats [2]:  number of garbage collections performed.  This can
00320  *        be excessively high if Alen is close to the minimum
00321  *        required value.
00322  *
00323  *    stats [3]:  status code.  < 0 is an error code.
00324  *          > 1 is a warning or notice.
00325  *
00326  *        0 OK.  Each column of the input matrix contained row
00327  *      indices in increasing order, with no duplicates.
00328  *
00329  *        1 OK, but columns of input matrix were jumbled (unsorted
00330  *      columns or duplicate entries).  CCOLAMD had to do some
00331  *      extra work to sort the matrix first and remove
00332  *      duplicate entries, but it still was able to return a
00333  *      valid permutation (return value of ccolamd was TRUE).
00334  *
00335  *      stats [4]: highest column index of jumbled columns
00336  *      stats [5]: last seen duplicate or unsorted row index
00337  *      stats [6]: number of duplicate or unsorted row indices
00338  *
00339  *        -1  A is a null pointer
00340  *
00341  *        -2  p is a null pointer
00342  *
00343  *        -3  n_row is negative.  stats [4]: n_row
00344  *
00345  *        -4  n_col is negative.  stats [4]: n_col
00346  *
00347  *        -5  number of nonzeros in matrix is negative
00348  *
00349  *      stats [4]: number of nonzeros, p [n_col]
00350  *
00351  *        -6  p [0] is nonzero
00352  *
00353  *      stats [4]: p [0]
00354  *
00355  *        -7  A is too small
00356  *
00357  *      stats [4]: required size
00358  *      stats [5]: actual size (Alen)
00359  *
00360  *        -8  a column has a negative number of entries
00361  *
00362  *      stats [4]: column with < 0 entries
00363  *      stats [5]: number of entries in col
00364  *
00365  *        -9  a row index is out of bounds
00366  *
00367  *      stats [4]: column with bad row index
00368  *      stats [5]: bad row index
00369  *      stats [6]: n_row, # of rows of matrx
00370  *
00371  *        -10 (unused; see csymamd)
00372  *
00373  *      int cmember [n_col] ;   Input argument.
00374  *
00375  *    cmember is new to CCOLAMD.  It did not appear in COLAMD.
00376  *    It places contraints on the output ordering.  s = cmember [j]
00377  *    gives the constraint set s that contains the column j
00378  *    (Restriction: 0 <= s < n_col).  In the output column
00379  *    permutation, all columns in set 0 appear first, followed by
00380  *    all columns in set 1, and so on.  If NULL, all columns are
00381  *    treated as if they were in a single constraint set, and you
00382  *    will obtain the same ordering as COLAMD (with one exception:
00383  *    the dense row/column threshold and other default knobs in
00384  *    CCOLAMD and COLAMD are different).
00385  *
00386  *  Example:
00387  *
00388  *      See
00389  *      http://www.cise.ufl.edu/research/sparse/ccolamd/ccolamd_example.c
00390  *      for a complete example.
00391  *
00392  *      To order the columns of a 5-by-4 matrix with 11 nonzero entries in
00393  *      the following nonzero pattern
00394  *
00395  *        x 0 x 0
00396  *    x 0 x x
00397  *    0 x x 0
00398  *    0 0 x x
00399  *    x x 0 0
00400  *
00401  *      with default knobs, no output statistics, and no ordering
00402  *      constraints, do the following:
00403  *
00404  *    #include "ccolamd.h"
00405  *    #define ALEN 144
00406  *    int A [ALEN] = {0, 1, 4, 2, 4, 0, 1, 2, 3, 1, 3} ;
00407  *    int p [ ] = {0, 3, 5, 9, 11} ;
00408  *    int stats [CCOLAMD_STATS] ;
00409  *    ccolamd (5, 4, ALEN, A, p, (double *) NULL, stats, NULL) ;
00410  *
00411  *      The permutation is returned in the array p, and A is destroyed.
00412  *
00413  *  ----------------------------------------------------------------------------
00414  *  csymamd:
00415  *  ----------------------------------------------------------------------------
00416  *
00417  *  C syntax:
00418  *
00419  *      #include "ccolamd.h"
00420  *
00421  *      int csymamd (int n, int *A, int *p, int *perm,
00422  *        double knobs [CCOLAMD_KNOBS], int stats [CCOLAMD_STATS],
00423  *    void (*allocate) (size_t, size_t), void (*release) (void *),
00424  *    int *cmember, int stype) ;
00425  *
00426  *      UF_long csymamd_l (UF_long n, UF_long *A, UF_long *p, UF_long *perm,
00427  *        double knobs [CCOLAMD_KNOBS], UF_long stats [CCOLAMD_STATS],
00428  *    void (*allocate) (size_t, size_t), void (*release) (void *),
00429  *    UF_long *cmember, UF_long stype) ;
00430  *
00431  *  Purpose:
00432  *
00433  *        The csymamd routine computes an ordering P of a symmetric sparse
00434  *      matrix A such that the Cholesky factorization PAP' = LL' remains
00435  *      sparse.  It is based on a column ordering of a matrix M constructed
00436  *      so that the nonzero pattern of M'M is the same as A.  Either the
00437  *      lower or upper triangular part of A can be used, or the pattern
00438  *      A+A' can be used.  You must pass your selected memory allocator
00439  *      (usually calloc/free or mxCalloc/mxFree) to csymamd, for it to
00440  *      allocate memory for the temporary matrix M.
00441  *
00442  *  Returns:
00443  *
00444  *      TRUE (1) if successful, FALSE (0) otherwise.
00445  *
00446  *  Arguments:
00447  *
00448  *      int n ;   Input argument.
00449  *
00450  *        Number of rows and columns in the symmetrix matrix A.
00451  *    Restriction:  n >= 0.
00452  *    csymamd returns FALSE if n is negative.
00453  *
00454  *      int A [nnz] ; Input argument.
00455  *
00456  *        A is an integer array of size nnz, where nnz = p [n].
00457  *
00458  *    The row indices of the entries in column c of the matrix are
00459  *    held in A [(p [c]) ... (p [c+1]-1)].  The row indices in a
00460  *    given column c need not be in ascending order, and duplicate
00461  *    row indices may be present.  However, csymamd will run faster
00462  *    if the columns are in sorted order with no duplicate entries.
00463  *
00464  *    The matrix is 0-based.  That is, rows are in the range 0 to
00465  *    n-1, and columns are in the range 0 to n-1.  csymamd
00466  *    returns FALSE if any row index is out of range.
00467  *
00468  *    The contents of A are not modified.
00469  *
00470  *      int p [n+1] ;     Input argument.
00471  *
00472  *    p is an integer array of size n+1.  On input, it holds the
00473  *    "pointers" for the column form of the matrix A.  Column c of
00474  *    the matrix A is held in A [(p [c]) ... (p [c+1]-1)].  The first
00475  *    entry, p [0], must be zero, and p [c] <= p [c+1] must hold
00476  *    for all c in the range 0 to n-1.  The value p [n] is
00477  *    thus the total number of entries in the pattern of the matrix A.
00478  *    csymamd returns FALSE if these conditions are not met.
00479  *
00480  *    The contents of p are not modified.
00481  *
00482  *      int perm [n+1] ;    Output argument.
00483  *
00484  *    On output, if csymamd returns TRUE, the array perm holds the
00485  *    permutation P, where perm [0] is the first index in the new
00486  *    ordering, and perm [n-1] is the last.  That is, perm [k] = j
00487  *    means that row and column j of A is the kth column in PAP',
00488  *    where k is in the range 0 to n-1 (perm [0] = j means
00489  *    that row and column j of A are the first row and column in
00490  *    PAP').  The array is used as a workspace during the ordering,
00491  *    which is why it must be of length n+1, not just n.
00492  *
00493  *      double knobs [CCOLAMD_KNOBS] ;  Input argument.
00494  *
00495  *    See colamd_set_defaults for a description.
00496  *
00497  *      int stats [CCOLAMD_STATS] ;   Output argument.
00498  *
00499  *    Statistics on the ordering, and error status.
00500  *    See ccolamd.h for related definitions.
00501  *    csymand returns FALSE if stats is not present.
00502  *
00503  *    stats [0]:  number of dense or empty row and columns ignored
00504  *        (and ordered last in the output permutation perm, subject
00505  *        to the constraints).  Note that a row/column can become
00506  *        "empty" if it contains only "dense" and/or "empty"
00507  *        columns/rows.
00508  *
00509  *    stats [1]:  (same as stats [0])
00510  *
00511  *    stats [2]:  number of garbage collections performed.
00512  *
00513  *    stats [3]:  status code.  < 0 is an error code.
00514  *          > 1 is a warning or notice.
00515  *
00516  *        0 to -9: same as ccolamd, with n replacing n_col and n_row,
00517  *      and -3 and -7 are unused.
00518  *
00519  *        -10 out of memory (unable to allocate temporary workspace
00520  *          for M or count arrays using the "allocate" routine
00521  *          passed into csymamd).
00522  *
00523  *      void * (*allocate) (size_t, size_t)
00524  *
00525  *        A pointer to a function providing memory allocation.  The
00526  *    allocated memory must be returned initialized to zero.  For a
00527  *    C application, this argument should normally be a pointer to
00528  *    calloc.  For a MATLAB mexFunction, the routine mxCalloc is
00529  *    passed instead.
00530  *
00531  *      void (*release) (size_t, size_t)
00532  *
00533  *        A pointer to a function that frees memory allocated by the
00534  *    memory allocation routine above.  For a C application, this
00535  *    argument should normally be a pointer to free.  For a MATLAB
00536  *    mexFunction, the routine mxFree is passed instead.
00537  *
00538  *      int cmember [n] ;   Input argument.
00539  *
00540  *    Same as ccolamd, except that cmember is of size n, and it places
00541  *    contraints symmetrically, on both the row and column ordering.
00542  *    Entries in cmember must be in the range 0 to n-1.
00543  *
00544  *      int stype ;     Input argument.
00545  *
00546  *    If stype < 0, then only the strictly lower triangular part of
00547  *    A is accessed.  The upper triangular part is assumed to be the
00548  *    transpose of the lower triangular part.  This is the same as
00549  *    SYMAMD, which did not have an stype parameter.
00550  *
00551  *    If stype > 0, only the strictly upper triangular part of A is
00552  *    accessed.  The lower triangular part is assumed to be the
00553  *    transpose of the upper triangular part.
00554  *
00555  *    If stype == 0, then the nonzero pattern of A+A' is ordered.
00556  *
00557  *  ----------------------------------------------------------------------------
00558  *  ccolamd_report:
00559  *  ----------------------------------------------------------------------------
00560  *
00561  *  C syntax:
00562  *
00563  *      #include "ccolamd.h"
00564  *      ccolamd_report (int stats [CCOLAMD_STATS]) ;
00565  *      ccolamd_l_report (UF_long stats [CCOLAMD_STATS]) ;
00566  *
00567  *  Purpose:
00568  *
00569  *      Prints the error status and statistics recorded in the stats
00570  *      array on the standard error output (for a standard C routine)
00571  *      or on the MATLAB output (for a mexFunction).
00572  *
00573  *  Arguments:
00574  *
00575  *      int stats [CCOLAMD_STATS] ; Input only.  Statistics from ccolamd.
00576  *
00577  *
00578  *  ----------------------------------------------------------------------------
00579  *  csymamd_report:
00580  *  ----------------------------------------------------------------------------
00581  *
00582  *  C syntax:
00583  *
00584  *      #include "ccolamd.h"
00585  *      csymamd_report (int stats [CCOLAMD_STATS]) ;
00586  *      csymamd_l_report (UF_long stats [CCOLAMD_STATS]) ;
00587  *
00588  *  Purpose:
00589  *
00590  *      Prints the error status and statistics recorded in the stats
00591  *      array on the standard error output (for a standard C routine)
00592  *      or on the MATLAB output (for a mexFunction).
00593  *
00594  *  Arguments:
00595  *
00596  *      int stats [CCOLAMD_STATS] ; Input only.  Statistics from csymamd.
00597  *
00598  */
00599
00600
00601 /* ========================================================================== */
00602 /* === Scaffolding code definitions  ======================================== */
00603 /* ========================================================================== */
00604
00605 /* Ensure that debugging is turned off: */
00606 #ifndef NDEBUG
00607 #define NDEBUG
00608 #endif
00609
00610 /* turn on debugging by uncommenting the following line
00611  #undef NDEBUG
00612  */
00613
00614 /* ========================================================================== */
00615 /* === Include files ======================================================== */
00616 /* ========================================================================== */
00617
00618 #include "amesos_ccolamd.h"
00619
00620 #include <stdlib.h>
00621 #include <math.h>
00622 #include <limits.h>
00623
00624 #ifdef MATLAB_MEX_FILE
00625 #include "mex.h"
00626 #include "matrix.h"
00627 #endif
00628
00629 #if !defined (NPRINT) || !defined (NDEBUG)
00630 #include <stdio.h>
00631 #endif
00632
00633 #ifndef NULL
00634 #define NULL ((void *) 0)
00635 #endif
00636
00637 /* ========================================================================== */
00638 /* === int or UF_long ======================================================= */
00639 /* ========================================================================== */
00640
00641 /* define UF_long */
00642 #include "amesos_UFconfig.h"
00643
00644 #ifdef DLONG
00645
00646 #define Int UF_long
00647 #define ID  UF_long_id
00648 #define Int_MAX UF_long_max
00649
00650 #define CCOLAMD_recommended amesos_ccolamd_l_recommended
00651 #define CCOLAMD_set_defaults amesos_ccolamd_l_set_defaults
00652 #define CCOLAMD_2 amesos_ccolamd2_l
00653 #define CCOLAMD_MAIN amesos_ccolamd_l
00654 #define CCOLAMD_apply_order amesos_ccolamd_l_apply_order
00655 #define CCOLAMD_postorder amesos_ccolamd_l_postorder
00656 #define CCOLAMD_post_tree amesos_ccolamd_l_post_tree
00657 #define CCOLAMD_fsize amesos_ccolamd_l_fsize
00658 #define CSYMAMD_MAIN amesos_csymamd_l
00659 #define CCOLAMD_report amesos_ccolamd_l_report
00660 #define CSYMAMD_report amesos_csymamd_l_report
00661
00662 #else
00663
00664 #define Int int
00665 #define ID "%d"
00666 #define Int_MAX INT_MAX
00667
00668 #define CCOLAMD_recommended amesos_ccolamd_recommended
00669 #define CCOLAMD_set_defaults amesos_ccolamd_set_defaults
00670 #define CCOLAMD_2 amesos_ccolamd2
00671 #define CCOLAMD_MAIN amesos_ccolamd
00672 #define CCOLAMD_apply_order amesos_ccolamd_apply_order
00673 #define CCOLAMD_postorder amesos_ccolamd_postorder
00674 #define CCOLAMD_post_tree amesos_ccolamd_post_tree
00675 #define CCOLAMD_fsize amesos_ccolamd_fsize
00676 #define CSYMAMD_MAIN amesos_csymamd
00677 #define CCOLAMD_report amesos_ccolamd_report
00678 #define CSYMAMD_report amesos_csymamd_report
00679
00680 #endif
00681
00682 /* ========================================================================== */
00683 /* === Row and Column structures ============================================ */
00684 /* ========================================================================== */
00685
00686 typedef struct CColamd_Col_struct
00687 {
00688     /* size of this struct is 8 integers if no padding occurs */
00689
00690     Int start ;   /* index for A of first row in this column, or DEAD */
00691       /* if column is dead */
00692     Int length ;  /* number of rows in this column */
00693     union
00694     {
00695   Int thickness ; /* number of original columns represented by this */
00696       /* col, if the column is alive */
00697   Int parent ;  /* parent in parent tree super-column structure, if */
00698       /* the column is dead */
00699     } shared1 ;
00700     union
00701     {
00702   Int score ;
00703   Int order ;
00704     } shared2 ;
00705     union
00706     {
00707   Int headhash ;  /* head of a hash bucket, if col is at the head of */
00708       /* a degree list */
00709   Int hash ;  /* hash value, if col is not in a degree list */
00710   Int prev ;  /* previous column in degree list, if col is in a */
00711       /* degree list (but not at the head of a degree list) */
00712     } shared3 ;
00713     union
00714     {
00715   Int degree_next ; /* next column, if col is in a degree list */
00716   Int hash_next ;   /* next column, if col is in a hash list */
00717     } shared4 ;
00718
00719     Int nextcol ;       /* next column in this supercolumn */
00720     Int lastcol ;       /* last column in this supercolumn */
00721
00722 } CColamd_Col ;
00723
00724
00725 typedef struct CColamd_Row_struct
00726 {
00727     /* size of this struct is 6 integers if no padding occurs */
00728
00729     Int start ;   /* index for A of first col in this row */
00730     Int length ;  /* number of principal columns in this row */
00731     union
00732     {
00733   Int degree ;  /* number of principal & non-principal columns in row */
00734   Int p ;   /* used as a row pointer in init_rows_cols () */
00735     } shared1 ;
00736     union
00737     {
00738   Int mark ;  /* for computing set differences and marking dead rows*/
00739   Int first_column ;/* first column in row (used in garbage collection) */
00740     } shared2 ;
00741
00742     Int thickness ;     /* number of original rows represented by this row */
00743                         /* that are not yet pivotal */
00744     Int front ;         /* -1 if an original row */
00745           /* k if this row represents the kth frontal matrix */
00746                         /* where k goes from 0 to at most n_col-1 */
00747
00748 } CColamd_Row ;
00749
00750 /* ========================================================================== */
00751 /* === basic definitions ==================================================== */
00752 /* ========================================================================== */
00753
00754 #define EMPTY (-1)
00755 #define MAX(a,b) (((a) > (b)) ? (a) : (b))
00756 #define MIN(a,b) (((a) < (b)) ? (a) : (b))
00757
00758 /* Routines are either PUBLIC (user-callable) or PRIVATE (not user-callable) */
00759 #define GLOBAL
00760 #define PUBLIC
00761 #define PRIVATE static
00762
00763 #define DENSE_DEGREE(alpha,n) \
00764     ((Int) MAX (16.0, (alpha) * sqrt ((double) (n))))
00765
00766 #define CMEMBER(c) ((cmember == (Int *) NULL) ? (0) : (cmember [c]))
00767
00768 /* True if x is NaN */
00769 #define SCALAR_IS_NAN(x)        ((x) != (x))
00770
00771 /* true if an integer (stored in double x) would overflow (or if x is NaN) */
00772 #define INT_OVERFLOW(x) ((!((x) * (1.0+1e-8) <= (double) Int_MAX)) \
00773                         || SCALAR_IS_NAN (x))
00774
00775 #define ONES_COMPLEMENT(r) (-(r)-1)
00776 #undef TRUE
00777 #undef FALSE
00778 #define TRUE (1)
00779 #define FALSE (0)
00780
00781 /* Row and column status */
00782 #define ALIVE (0)
00783 #define DEAD  (-1)
00784
00785 /* Column status */
00786 #define DEAD_PRINCIPAL    (-1)
00787 #define DEAD_NON_PRINCIPAL  (-2)
00788
00789 /* Macros for row and column status update and checking. */
00790 #define ROW_IS_DEAD(r)      ROW_IS_MARKED_DEAD (Row[r].shared2.mark)
00791 #define ROW_IS_MARKED_DEAD(row_mark)  (row_mark < ALIVE)
00792 #define ROW_IS_ALIVE(r)     (Row [r].shared2.mark >= ALIVE)
00793 #define COL_IS_DEAD(c)      (Col [c].start < ALIVE)
00794 #define COL_IS_ALIVE(c)     (Col [c].start >= ALIVE)
00795 #define COL_IS_DEAD_PRINCIPAL(c)  (Col [c].start == DEAD_PRINCIPAL)
00796 #define KILL_ROW(r)     { Row [r].shared2.mark = DEAD ; }
00797 #define KILL_PRINCIPAL_COL(c)   { Col [c].start = DEAD_PRINCIPAL ; }
00798 #define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; }
00799
00800
00801 /* ========================================================================== */
00802 /* === ccolamd reporting mechanism ========================================== */
00803 /* ========================================================================== */
00804
00805 #if defined (MATLAB_MEX_FILE) || defined (MATHWORKS)
00806 /* In MATLAB, matrices are 1-based to the user, but 0-based internally */
00807 #define INDEX(i) ((i)+1)
00808 #else
00809 /* In C, matrices are 0-based and indices are reported as such in *_report */
00810 #define INDEX(i) (i)
00811 #endif
00812
00813 /* All output goes through the PRINTF macro.  */
00814 #define PRINTF(params) { if (amesos_ccolamd_printf != NULL) (void) amesos_ccolamd_printf params ; }
00815
00816
00817 /* ========================================================================== */
00818 /* === Debugging prototypes and definitions ================================= */
00819 /* ========================================================================== */
00820
00821 #ifndef NDEBUG
00822
00823 #include <assert.h>
00824
00825 /* debug print level, present only when debugging */
00826 PRIVATE Int ccolamd_debug ;
00827
00828 /* debug print statements */
00829 #define DEBUG0(params) { PRINTF (params) ; }
00830 #define DEBUG1(params) { if (ccolamd_debug >= 1) PRINTF (params) ; }
00831 #define DEBUG2(params) { if (ccolamd_debug >= 2) PRINTF (params) ; }
00832 #define DEBUG3(params) { if (ccolamd_debug >= 3) PRINTF (params) ; }
00833 #define DEBUG4(params) { if (ccolamd_debug >= 4) PRINTF (params) ; }
00834
00835 #ifdef MATLAB_MEX_FILE
00836 #define ASSERT(expression) (mxAssert ((expression), ""))
00837 #else
00838 #define ASSERT(expression) (assert (expression))
00839 #endif
00840
00841 PRIVATE void ccolamd_get_debug
00842 (
00843     char *method
00844 ) ;
00845
00846 PRIVATE void debug_mark
00847 (
00848     Int n_row,
00849     CColamd_Row Row [],
00850     Int tag_mark,
00851     Int max_mark
00852 ) ;
00853
00854 PRIVATE void debug_matrix
00855 (
00856     Int n_row,
00857     Int n_col,
00858     CColamd_Row Row [],
00859     CColamd_Col Col [],
00860     Int A []
00861 ) ;
00862
00863 PRIVATE void debug_structures
00864 (
00865     Int n_row,
00866     Int n_col,
00867     CColamd_Row Row [],
00868     CColamd_Col Col [],
00869     Int A [],
00870     Int in_cset [],
00871     Int cset_start []
00872 ) ;
00873
00874 PRIVATE void dump_super
00875 (
00876     Int super_c,
00877     CColamd_Col Col [],
00878     Int n_col
00879 ) ;
00880
00881 PRIVATE void debug_deg_lists
00882 (
00883     Int n_row,
00884     Int n_col,
00885     CColamd_Row Row [ ],
00886     CColamd_Col Col [ ],
00887     Int head [ ],
00888     Int min_score,
00889     Int should,
00890     Int max_deg
00891 ) ;
00892
00893 #else
00894
00895 /* === No debugging ========================================================= */
00896
00897 #define DEBUG0(params) ;
00898 #define DEBUG1(params) ;
00899 #define DEBUG2(params) ;
00900 #define DEBUG3(params) ;
00901 #define DEBUG4(params) ;
00902
00903 #define ASSERT(expression)
00904
00905 #endif
00906
00907 /* ========================================================================== */
00908 /* === Prototypes of PRIVATE routines ======================================= */
00909 /* ========================================================================== */
00910
00911 PRIVATE Int init_rows_cols
00912 (
00913     Int n_row,
00914     Int n_col,
00915     CColamd_Row Row [ ],
00916     CColamd_Col Col [ ],
00917     Int A [ ],
00918     Int p [ ],
00919     Int stats [CCOLAMD_STATS]
00920 ) ;
00921
00922 PRIVATE void init_scoring
00923 (
00924     Int n_row,
00925     Int n_col,
00926     CColamd_Row Row [ ],
00927     CColamd_Col Col [ ],
00928     Int A [ ],
00929     Int head [ ],
00930     double knobs [CCOLAMD_KNOBS],
00931     Int *p_n_row2,
00932     Int *p_n_col2,
00933     Int *p_max_deg,
00934     Int cmember [ ],
00935     Int n_cset,
00936     Int cset_start [ ],
00937     Int dead_cols [ ],
00938     Int *p_ndense_row,    /* number of dense rows */
00939     Int *p_nempty_row,    /* number of original empty rows */
00940     Int *p_nnewlyempty_row, /* number of newly empty rows */
00941     Int *p_ndense_col,    /* number of dense cols (excl "empty" cols) */
00942     Int *p_nempty_col,    /* number of original empty cols */
00943     Int *p_nnewlyempty_col  /* number of newly empty cols */
00944 ) ;
00945
00946 PRIVATE Int find_ordering
00947 (
00948     Int n_row,
00949     Int n_col,
00950     Int Alen,
00951     CColamd_Row Row [ ],
00952     CColamd_Col Col [ ],
00953     Int A [ ],
00954     Int head [ ],
00955 #ifndef NDEBUG
00956     Int n_col2,
00957 #endif
00958     Int max_deg,
00959     Int pfree,
00960     Int cset [ ],
00961     Int cset_start [ ],
00962 #ifndef NDEBUG
00963     Int n_cset,
00964 #endif
00965     Int cmember [ ],
00966     Int Front_npivcol [ ],
00967     Int Front_nrows [ ],
00968     Int Front_ncols [ ],
00969     Int Front_parent [ ],
00970     Int Front_cols [ ],
00971     Int *p_nfr,
00972     Int aggressive,
00973     Int InFront [ ],
00974     Int order_for_lu
00975 ) ;
00976
00977 PRIVATE void detect_super_cols
00978 (
00979 #ifndef NDEBUG
00980     Int n_col,
00981     CColamd_Row Row [ ],
00982 #endif
00983     CColamd_Col Col [ ],
00984     Int A [ ],
00985     Int head [ ],
00986     Int row_start,
00987     Int row_length,
00988     Int in_set [ ]
00989 ) ;
00990
00991 PRIVATE Int garbage_collection
00992 (
00993     Int n_row,
00994     Int n_col,
00995     CColamd_Row Row [ ],
00996     CColamd_Col Col [ ],
00997     Int A [ ],
00998     Int *pfree
00999 ) ;
01000
01001 PRIVATE Int clear_mark
01002 (
01003     Int tag_mark,
01004     Int max_mark,
01005     Int n_row,
01006     CColamd_Row Row [ ]
01007 ) ;
01008
01009 PRIVATE void print_report
01010 (
01011     char *method,
01012     Int stats [CCOLAMD_STATS]
01013 ) ;
01014
01015
01016 /* ========================================================================== */
01017 /* === USER-CALLABLE ROUTINES: ============================================== */
01018 /* ========================================================================== */
01019
01020
01021 /* ========================================================================== */
01022 /* === ccolamd_recommended ================================================== */
01023 /* ========================================================================== */
01024
01025 /*
01026  *  The ccolamd_recommended routine returns the suggested size for Alen.  This
01027  *  value has been determined to provide good balance between the number of
01028  *  garbage collections and the memory requirements for ccolamd.  If any
01029  *  argument is negative, or if integer overflow occurs, a 0 is returned as
01030  *  an error condition.
01031  *
01032  *  2*nnz space is required for the row and column indices of the matrix
01033  *  (or 4*n_col, which ever is larger).
01034  *
01035  *  CCOLAMD_C (n_col) + CCOLAMD_R (n_row) space is required for the Col and Row
01036  *  arrays, respectively, which are internal to ccolamd.  This is equal to
01037  *  8*n_col + 6*n_row if the structures are not padded.
01038  *
01039  *  An additional n_col space is the minimal amount of "elbow room",
01040  *  and nnz/5 more space is recommended for run time efficiency.
01041  *
01042  *  The remaining (((3 * n_col) + 1) + 5 * (n_col + 1) + n_row) space is
01043  *  for other workspace used in ccolamd which did not appear in colamd.
01044  */
01045
01046 /* add two values of type size_t, and check for integer overflow */
01047 static size_t t_add (size_t a, size_t b, int *ok)
01048 {
01049     (*ok) = (*ok) && ((a + b) >= MAX (a,b)) ;
01050     return ((*ok) ? (a + b) : 0) ;
01051 }
01052
01053 /* compute a*k where k is a small integer, and check for integer overflow */
01054 static size_t t_mult (size_t a, size_t k, int *ok)
01055 {
01056     size_t i, s = 0 ;
01057     for (i = 0 ; i < k ; i++)
01058     {
01059   s = t_add (s, a, ok) ;
01060     }
01061     return (s) ;
01062 }
01063
01064 /* size of the Col and Row structures */
01065 #define CCOLAMD_C(n_col,ok) \
01066     ((t_mult (t_add (n_col, 1, ok), sizeof (CColamd_Col), ok) / sizeof (Int)))
01067
01068 #define CCOLAMD_R(n_row,ok) \
01069     ((t_mult (t_add (n_row, 1, ok), sizeof (CColamd_Row), ok) / sizeof (Int)))
01070
01071 /*
01072 #define CCOLAMD_RECOMMENDED(nnz, n_row, n_col) \
01073       MAX (2 * nnz, 4 * n_col) + \
01074       CCOLAMD_C (n_col) + CCOLAMD_R (n_row) + n_col + (nnz / 5) \
01075       + ((3 * n_col) + 1) + 5 * (n_col + 1) + n_row
01076  */
01077
01078 static size_t ccolamd_need (Int nnz, Int n_row, Int n_col, int *ok)
01079 {
01080
01081     /* ccolamd_need, compute the following, and check for integer overflow:
01082   need = MAX (2*nnz, 4*n_col) + n_col +
01083     Col_size + Row_size +
01084     (3*n_col+1) + (5*(n_col+1)) + n_row ;
01085     */
01086     size_t s, c, r, t ;
01087
01088     /* MAX (2*nnz, 4*n_col) */
01089     s = t_mult (nnz, 2, ok) ;     /* 2*nnz */
01090     t = t_mult (n_col, 4, ok) ;     /* 4*n_col */
01091     s = MAX (s,t) ;
01092
01093     s = t_add (s, n_col, ok) ;      /* bare minimum elbow room */
01094
01095     /* Col and Row arrays */
01096     c = CCOLAMD_C (n_col, ok) ;     /* size of column structures */
01097     r = CCOLAMD_R (n_row, ok) ;     /* size of row structures */
01098     s = t_add (s, c, ok) ;
01099     s = t_add (s, r, ok) ;
01100
01101     c = t_mult (n_col, 3, ok) ;     /* 3*n_col + 1 */
01102     c = t_add (c, 1, ok) ;
01103     s = t_add (s, c, ok) ;
01104
01105     c = t_add (n_col, 1, ok) ;      /* 5 * (n_col + 1) */
01106     c = t_mult (c, 5, ok) ;
01107     s = t_add (s, c, ok) ;
01108
01109     s = t_add (s, n_row, ok) ;      /* n_row */
01110
01111     return (ok ? s : 0) ;
01112 }
01113
01114 PUBLIC size_t CCOLAMD_recommended /* returns recommended value of Alen. */
01115 (
01116     /* === Parameters ======================================================= */
01117
01118     Int nnz,      /* number of nonzeros in A */
01119     Int n_row,      /* number of rows in A */
01120     Int n_col     /* number of columns in A */
01121 )
01122 {
01123     size_t s ;
01124     int ok = TRUE ;
01125     if (nnz < 0 || n_row < 0 || n_col < 0)
01126     {
01127   return (0) ;
01128     }
01129     s = ccolamd_need (nnz, n_row, n_col, &ok) ; /* bare minimum needed */
01130     s = t_add (s, nnz/5, &ok) ;     /* extra elbow room */
01131     ok = ok && (s < Int_MAX) ;
01132     return (ok ? s : 0) ;
01133 }
01134
01135
01136 /* ========================================================================== */
01137 /* === ccolamd_set_defaults ================================================= */
01138 /* ========================================================================== */
01139
01140 /*
01141  *  The ccolamd_set_defaults routine sets the default values of the user-
01142  *  controllable parameters for ccolamd.
01143  */
01144
01145 PUBLIC void CCOLAMD_set_defaults
01146 (
01147     /* === Parameters ======================================================= */
01148
01149     double knobs [CCOLAMD_KNOBS]    /* knob array */
01150 )
01151 {
01152     /* === Local variables ================================================== */
01153
01154     Int i ;
01155
01156     if (!knobs)
01157     {
01158   return ;      /* no knobs to initialize */
01159     }
01160     for (i = 0 ; i < CCOLAMD_KNOBS ; i++)
01161     {
01162   knobs [i] = 0 ;
01163     }
01164     knobs [CCOLAMD_DENSE_ROW] = 10 ;
01165     knobs [CCOLAMD_DENSE_COL] = 10 ;
01166     knobs [CCOLAMD_AGGRESSIVE] = TRUE ; /* default: do aggressive absorption*/
01167     knobs [CCOLAMD_LU] = FALSE ;  /* default: order for Cholesky */
01168 }
01169
01170
01171 /* ========================================================================== */
01172 /* === symamd =============================================================== */
01173 /* ========================================================================== */
01174
01175 PUBLIC Int CSYMAMD_MAIN   /* return TRUE if OK, FALSE otherwise */
01176 (
01177     /* === Parameters ======================================================= */
01178
01179     Int n,        /* number of rows and columns of A */
01180     Int A [ ],        /* row indices of A */
01181     Int p [ ],        /* column pointers of A */
01182     Int perm [ ],     /* output permutation, size n+1 */
01183     double knobs [CCOLAMD_KNOBS], /* parameters (uses defaults if NULL) */
01184     Int stats [CCOLAMD_STATS],    /* output statistics and error codes */
01185     void * (*allocate) (size_t, size_t),/* pointer to calloc (ANSI C) or */
01186           /* mxCalloc (for MATLAB mexFunction) */
01187     void (*release) (void *),   /* pointer to free (ANSI C) or */
01188               /* mxFree (for MATLAB mexFunction) */
01189     Int cmember [ ],      /* constraint set */
01190     Int stype             /* stype of A */
01191 )
01192 {
01193     /* === Local variables ================================================== */
01194
01195     double cknobs [CCOLAMD_KNOBS] ;
01196     double default_knobs [CCOLAMD_KNOBS] ;
01197
01198     Int *count ;    /* length of each column of M, and col pointer*/
01199     Int *mark ;     /* mark array for finding duplicate entries */
01200     Int *M ;      /* row indices of matrix M */
01201     size_t Mlen ;   /* length of M */
01202     Int n_row ;     /* number of rows in M */
01203     Int nnz ;     /* number of entries in A */
01204     Int i ;     /* row index of A */
01205     Int j ;     /* column index of A */
01206     Int k ;     /* row index of M */
01207     Int mnz ;     /* number of nonzeros in M */
01208     Int pp ;      /* index into a column of A */
01209     Int last_row ;    /* last row seen in the current column */
01210     Int length ;    /* number of nonzeros in a column */
01211     Int both ;      /* TRUE if ordering A+A' */
01212     Int upper ;     /* TRUE if ordering triu(A)+triu(A)' */
01213     Int lower ;     /* TRUE if ordering tril(A)+tril(A)' */
01214
01215 #ifndef NDEBUG
01216     ccolamd_get_debug ("csymamd") ;
01217 #endif
01218
01219     both = (stype == 0) ;
01220     upper = (stype > 0) ;
01221     lower = (stype < 0) ;
01222
01223     /* === Check the input arguments ======================================== */
01224
01225     if (!stats)
01226     {
01227   DEBUG1 (("csymamd: stats not present\n")) ;
01228   return (FALSE) ;
01229     }
01230     for (i = 0 ; i < CCOLAMD_STATS ; i++)
01231     {
01232   stats [i] = 0 ;
01233     }
01234     stats [CCOLAMD_STATUS] = CCOLAMD_OK ;
01235     stats [CCOLAMD_INFO1] = -1 ;
01236     stats [CCOLAMD_INFO2] = -1 ;
01237
01238     if (!A)
01239     {
01240       stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_not_present ;
01241   DEBUG1 (("csymamd: A not present\n")) ;
01242   return (FALSE) ;
01243     }
01244
01245     if (!p)   /* p is not present */
01246     {
01247   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p_not_present ;
01248   DEBUG1 (("csymamd: p not present\n")) ;
01249       return (FALSE) ;
01250     }
01251
01252     if (n < 0)    /* n must be >= 0 */
01253     {
01254   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_ncol_negative ;
01255   stats [CCOLAMD_INFO1] = n ;
01256   DEBUG1 (("csymamd: n negative "ID" \n", n)) ;
01257       return (FALSE) ;
01258     }
01259
01260     nnz = p [n] ;
01261     if (nnz < 0)  /* nnz must be >= 0 */
01262     {
01263   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nnz_negative ;
01264   stats [CCOLAMD_INFO1] = nnz ;
01265   DEBUG1 (("csymamd: number of entries negative "ID" \n", nnz)) ;
01266   return (FALSE) ;
01267     }
01268
01269     if (p [0] != 0)
01270     {
01271   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p0_nonzero ;
01272   stats [CCOLAMD_INFO1] = p [0] ;
01273   DEBUG1 (("csymamd: p[0] not zero "ID"\n", p [0])) ;
01274   return (FALSE) ;
01275     }
01276
01277     /* === If no knobs, set default knobs =================================== */
01278
01279     if (!knobs)
01280     {
01281   CCOLAMD_set_defaults (default_knobs) ;
01282   knobs = default_knobs ;
01283     }
01284
01285     /* === Allocate count and mark ========================================== */
01286
01287     count = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01288     if (!count)
01289     {
01290   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
01291   DEBUG1 (("csymamd: allocate count (size "ID") failed\n", n+1)) ;
01292   return (FALSE) ;
01293     }
01294
01295     mark = (Int *) ((*allocate) (n+1, sizeof (Int))) ;
01296     if (!mark)
01297     {
01298   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
01299   (*release) ((void *) count) ;
01300   DEBUG1 (("csymamd: allocate mark (size "ID") failed\n", n+1)) ;
01301   return (FALSE) ;
01302     }
01303
01304     /* === Compute column counts of M, check if A is valid ================== */
01305
01306     stats [CCOLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
01307
01308     for (i = 0 ; i < n ; i++)
01309     {
01310       mark [i] = -1 ;
01311     }
01312     for (j = 0 ; j < n ; j++)
01313     {
01314   last_row = -1 ;
01315
01316   length = p [j+1] - p [j] ;
01317   if (length < 0)
01318   {
01319       /* column pointers must be non-decreasing */
01320       stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_col_length_negative ;
01321       stats [CCOLAMD_INFO1] = j ;
01322       stats [CCOLAMD_INFO2] = length ;
01323       (*release) ((void *) count) ;
01324       (*release) ((void *) mark) ;
01325       DEBUG1 (("csymamd: col "ID" negative length "ID"\n", j, length)) ;
01326       return (FALSE) ;
01327   }
01328
01329   for (pp = p [j] ; pp < p [j+1] ; pp++)
01330   {
01331       i = A [pp] ;
01332       if (i < 0 || i >= n)
01333       {
01334     /* row index i, in column j, is out of bounds */
01335     stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_row_index_out_of_bounds ;
01336     stats [CCOLAMD_INFO1] = j ;
01337     stats [CCOLAMD_INFO2] = i ;
01338     stats [CCOLAMD_INFO3] = n ;
01339     (*release) ((void *) count) ;
01340     (*release) ((void *) mark) ;
01341     DEBUG1 (("csymamd: row "ID" col "ID" out of bounds\n", i, j)) ;
01342     return (FALSE) ;
01343       }
01344
01345       if (i <= last_row || mark [i] == j)
01346       {
01347     /* row index is unsorted or repeated (or both), thus col */
01348     /* is jumbled.  This is a notice, not an error condition. */
01349     stats [CCOLAMD_STATUS] = CCOLAMD_OK_BUT_JUMBLED ;
01350     stats [CCOLAMD_INFO1] = j ;
01351     stats [CCOLAMD_INFO2] = i ;
01352     (stats [CCOLAMD_INFO3]) ++ ;
01353     DEBUG1 (("csymamd: row "ID" col "ID" unsorted/dupl.\n", i, j)) ;
01354       }
01355
01356       if (mark [i] != j)
01357       {
01358     if ((both && i != j) || (lower && i > j) || (upper && i < j))
01359     {
01360         /* row k of M will contain column indices i and j */
01361         count [i]++ ;
01362         count [j]++ ;
01363     }
01364       }
01365
01366       /* mark the row as having been seen in this column */
01367       mark [i] = j ;
01368
01369       last_row = i ;
01370   }
01371     }
01372
01373     /* === Compute column pointers of M ===================================== */
01374
01375     /* use output permutation, perm, for column pointers of M */
01376     perm [0] = 0 ;
01377     for (j = 1 ; j <= n ; j++)
01378     {
01379   perm [j] = perm [j-1] + count [j-1] ;
01380     }
01381     for (j = 0 ; j < n ; j++)
01382     {
01383   count [j] = perm [j] ;
01384     }
01385
01386     /* === Construct M ====================================================== */
01387
01388     mnz = perm [n] ;
01389     n_row = mnz / 2 ;
01390     Mlen = CCOLAMD_recommended (mnz, n_row, n) ;
01391     M = (Int *) ((*allocate) (Mlen, sizeof (Int))) ;
01392     DEBUG1 (("csymamd: M is "ID"-by-"ID" with "ID" entries, Mlen = %g\n",
01393       n_row, n, mnz, (double) Mlen)) ;
01394
01395     if (!M)
01396     {
01397   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_out_of_memory ;
01398   (*release) ((void *) count) ;
01399   (*release) ((void *) mark) ;
01400   DEBUG1 (("csymamd: allocate M (size %g) failed\n", (double) Mlen)) ;
01401   return (FALSE) ;
01402     }
01403
01404     k = 0 ;
01405
01406     if (stats [CCOLAMD_STATUS] == CCOLAMD_OK)
01407     {
01408   /* Matrix is OK */
01409   for (j = 0 ; j < n ; j++)
01410   {
01411       ASSERT (p [j+1] - p [j] >= 0) ;
01412       for (pp = p [j] ; pp < p [j+1] ; pp++)
01413       {
01414     i = A [pp] ;
01415     ASSERT (i >= 0 && i < n) ;
01416     if ((both && i != j) || (lower && i > j) || (upper && i < j))
01417     {
01418         /* row k of M contains column indices i and j */
01419         M [count [i]++] = k ;
01420         M [count [j]++] = k ;
01421         k++ ;
01422     }
01423       }
01424   }
01425     }
01426     else
01427     {
01428   /* Matrix is jumbled.  Do not add duplicates to M.  Unsorted cols OK. */
01429   DEBUG1 (("csymamd: Duplicates in A.\n")) ;
01430   for (i = 0 ; i < n ; i++)
01431   {
01432       mark [i] = -1 ;
01433   }
01434   for (j = 0 ; j < n ; j++)
01435   {
01436       ASSERT (p [j+1] - p [j] >= 0) ;
01437       for (pp = p [j] ; pp < p [j+1] ; pp++)
01438       {
01439     i = A [pp] ;
01440     ASSERT (i >= 0 && i < n) ;
01441     if (mark [i] != j)
01442     {
01443         if ((both && i != j) || (lower && i > j) || (upper && i<j))
01444         {
01445       /* row k of M contains column indices i and j */
01446       M [count [i]++] = k ;
01447       M [count [j]++] = k ;
01448       k++ ;
01449       mark [i] = j ;
01450         }
01451     }
01452       }
01453   }
01454     }
01455
01456     /* count and mark no longer needed */
01457     (*release) ((void *) mark) ;
01458     (*release) ((void *) count) ;
01459     ASSERT (k == n_row) ;
01460
01461     /* === Adjust the knobs for M =========================================== */
01462
01463     for (i = 0 ; i < CCOLAMD_KNOBS ; i++)
01464     {
01465   cknobs [i] = knobs [i] ;
01466     }
01467
01468     /* there are no dense rows in M */
01469     cknobs [CCOLAMD_DENSE_ROW] = -1 ;
01470     cknobs [CCOLAMD_DENSE_COL] = knobs [CCOLAMD_DENSE_ROW] ;
01471
01472     /* ensure CCSYMAMD orders for Cholesky, not LU */
01473     cknobs [CCOLAMD_LU] = FALSE ;
01474
01475     /* === Order the columns of M =========================================== */
01476
01477     (void) CCOLAMD_2 (n_row, n, (Int) Mlen, M, perm, cknobs, stats,
01478              (Int *) NULL, (Int *) NULL, (Int *) NULL, (Int *) NULL,
01479              (Int *) NULL, (Int *) NULL, (Int *) NULL, cmember) ;
01480
01481     /* === adjust statistics ================================================ */
01482
01483     /* a dense column in ccolamd means a dense row and col in csymamd */
01484     stats [CCOLAMD_DENSE_ROW] = stats [CCOLAMD_DENSE_COL] ;
01485
01486     /* === Free M =========================================================== */
01487
01488     (*release) ((void *) M) ;
01489     DEBUG1 (("csymamd: done.\n")) ;
01490     return (TRUE) ;
01491 }
01492
01493
01494 /* ========================================================================== */
01495 /* === ccolamd ============================================================== */
01496 /* ========================================================================== */
01497
01498 /*
01499  *  The colamd routine computes a column ordering Q of a sparse matrix
01500  *  A such that the LU factorization P(AQ) = LU remains sparse, where P is
01501  *  selected via partial pivoting.   The routine can also be viewed as
01502  *  providing a permutation Q such that the Cholesky factorization
01503  *  (AQ)'(AQ) = LL' remains sparse.
01504  */
01505
01506 PUBLIC Int CCOLAMD_MAIN
01507 (
01508     /* === Parameters ======================================================= */
01509
01510     Int n_row,      /* number of rows in A */
01511     Int n_col,      /* number of columns in A */
01512     Int Alen,     /* length of A */
01513     Int A [ ],      /* row indices of A */
01514     Int p [ ],      /* pointers to columns in A */
01515     double knobs [CCOLAMD_KNOBS],/* parameters (uses defaults if NULL) */
01516     Int stats [CCOLAMD_STATS],  /* output statistics and error codes */
01517     Int cmember [ ]   /* constraint set of A */
01518 )
01519 {
01520      return (CCOLAMD_2 (n_row, n_col, Alen, A, p, knobs, stats,
01521              (Int *) NULL, (Int *) NULL, (Int *) NULL, (Int *) NULL,
01522              (Int *) NULL, (Int *) NULL, (Int *) NULL, cmember)) ;
01523 }
01524
01525
01526 /* ========================================================================== */
01527 /* === ccolamd2 ============================================================= */
01528 /* ========================================================================== */
01529
01530 /* Identical to ccolamd, except that additional information about each frontal
01531  * matrix is returned to the caller.  Not intended to be directly called by
01532  * the user.
01533  */
01534
01535 PUBLIC Int CCOLAMD_2      /* returns TRUE if successful, FALSE otherwise */
01536 (
01537     /* === Parameters ======================================================= */
01538
01539     Int n_row,      /* number of rows in A */
01540     Int n_col,      /* number of columns in A */
01541     Int Alen,     /* length of A */
01542     Int A [ ],      /* row indices of A */
01543     Int p [ ],      /* pointers to columns in A */
01544     double knobs [CCOLAMD_KNOBS],/* parameters (uses defaults if NULL) */
01545     Int stats [CCOLAMD_STATS],  /* output statistics and error codes */
01546
01547     /* each Front array is of size n_col+1. */
01548     Int Front_npivcol [ ],  /* # pivot cols in each front */
01549     Int Front_nrows [ ],  /* # of rows in each front (incl. pivot rows) */
01550     Int Front_ncols [ ],  /* # of cols in each front (incl. pivot cols) */
01551     Int Front_parent [ ], /* parent of each front */
01552     Int Front_cols [ ],   /* link list of pivot columns for each front */
01553     Int *p_nfr,     /* total number of frontal matrices */
01554     Int InFront [ ],    /* InFront [row] = f if the original row was
01555          * absorbed into front f.  EMPTY if the row was
01556          * empty, dense, or not absorbed.  This array
01557          * has size n_row+1 */
01558     Int cmember [ ]   /* constraint set of A */
01559 )
01560 {
01561     /* === Local variables ================================================== */
01562
01563     Int i ;     /* loop index */
01564     Int nnz ;     /* nonzeros in A */
01565     size_t Row_size ;   /* size of Row [ ], in integers */
01566     size_t Col_size ;   /* size of Col [ ], in integers */
01567     size_t need ;   /* minimum required length of A */
01568     CColamd_Row *Row ;    /* pointer into A of Row [0..n_row] array */
01569     CColamd_Col *Col ;    /* pointer into A of Col [0..n_col] array */
01570     Int n_col2 ;    /* number of non-dense, non-empty columns */
01571     Int n_row2 ;    /* number of non-dense, non-empty rows */
01572     Int ngarbage ;    /* number of garbage collections performed */
01573     Int max_deg ;   /* maximum row degree */
01574     double default_knobs [CCOLAMD_KNOBS] ;  /* default knobs array */
01575
01576     Int n_cset ;    /* number of constraint sets */
01577     Int *cset ;     /* cset of A */
01578     Int *cset_start ;   /* pointer into cset */
01579     Int *temp_cstart ;    /* temp pointer to start of cset */
01580     Int *csize ;    /* temp pointer to cset size */
01581     Int ap ;      /* column index */
01582     Int order_for_lu ;    /* TRUE: order for LU, FALSE: for Cholesky */
01583
01584     Int ndense_row, nempty_row, parent, ndense_col,
01585       nempty_col, k, col, nfr, *Front_child, *Front_sibling, *Front_stack,
01586       *Front_order, *Front_size ;
01587     Int nnewlyempty_col, nnewlyempty_row ;
01588     Int aggressive ;
01589     Int row ;
01590     Int *dead_cols ;
01591     Int set1 ;
01592     Int set2 ;
01593     Int cs ;
01594
01595     int ok ;
01596
01597 #ifndef NDEBUG
01598     ccolamd_get_debug ("ccolamd") ;
01599 #endif
01600
01601     /* === Check the input arguments ======================================== */
01602
01603     if (!stats)
01604     {
01605   DEBUG1 (("ccolamd: stats not present\n")) ;
01606   return (FALSE) ;
01607     }
01608     for (i = 0 ; i < CCOLAMD_STATS ; i++)
01609     {
01610   stats [i] = 0 ;
01611     }
01612     stats [CCOLAMD_STATUS] = CCOLAMD_OK ;
01613     stats [CCOLAMD_INFO1] = -1 ;
01614     stats [CCOLAMD_INFO2] = -1 ;
01615
01616     if (!A)   /* A is not present */
01617     {
01618   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_not_present ;
01619   DEBUG1 (("ccolamd: A not present\n")) ;
01620   return (FALSE) ;
01621     }
01622
01623     if (!p)   /* p is not present */
01624     {
01625   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p_not_present ;
01626   DEBUG1 (("ccolamd: p not present\n")) ;
01627       return (FALSE) ;
01628     }
01629
01630     if (n_row < 0)  /* n_row must be >= 0 */
01631     {
01632   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nrow_negative ;
01633   stats [CCOLAMD_INFO1] = n_row ;
01634   DEBUG1 (("ccolamd: nrow negative "ID"\n", n_row)) ;
01635       return (FALSE) ;
01636     }
01637
01638     if (n_col < 0)  /* n_col must be >= 0 */
01639     {
01640   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_ncol_negative ;
01641   stats [CCOLAMD_INFO1] = n_col ;
01642   DEBUG1 (("ccolamd: ncol negative "ID"\n", n_col)) ;
01643       return (FALSE) ;
01644     }
01645
01646     nnz = p [n_col] ;
01647     if (nnz < 0)  /* nnz must be >= 0 */
01648     {
01649   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_nnz_negative ;
01650   stats [CCOLAMD_INFO1] = nnz ;
01651   DEBUG1 (("ccolamd: number of entries negative "ID"\n", nnz)) ;
01652   return (FALSE) ;
01653     }
01654
01655     if (p [0] != 0)
01656     {
01657   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_p0_nonzero ;
01658   stats [CCOLAMD_INFO1] = p [0] ;
01659   DEBUG1 (("ccolamd: p[0] not zero "ID"\n", p [0])) ;
01660   return (FALSE) ;
01661     }
01662
01663     /* === If no knobs, set default knobs =================================== */
01664
01665     if (!knobs)
01666     {
01667   CCOLAMD_set_defaults (default_knobs) ;
01668   knobs = default_knobs ;
01669     }
01670
01671     aggressive = (knobs [CCOLAMD_AGGRESSIVE] != FALSE) ;
01672     order_for_lu = (knobs [CCOLAMD_LU] != FALSE) ;
01673
01674     /* === Allocate workspace from array A ================================== */
01675
01676     ok = TRUE ;
01677     Col_size = CCOLAMD_C (n_col, &ok) ;
01678     Row_size = CCOLAMD_R (n_row, &ok) ;
01679
01680     /* min size of A is 2nnz+ncol.  cset and cset_start are of size 2ncol+1 */
01681     /* Each of the 5 fronts is of size n_col + 1. InFront is of size nrow.  */
01682
01683     /*
01684     need = MAX (2*nnz, 4*n_col) + n_col +
01685         Col_size + Row_size +
01686     (3*n_col+1) + (5*(n_col+1)) + n_row ;
01687     */
01688     need = ccolamd_need (nnz, n_row, n_col, &ok) ;
01689
01690     if (!ok || need > (size_t) Alen || need > Int_MAX)
01691     {
01692   /* not enough space in array A to perform the ordering */
01693   stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_A_too_small ;
01694   stats [CCOLAMD_INFO1] = need ;
01695   stats [CCOLAMD_INFO2] = Alen ;
01696   DEBUG1 (("ccolamd: Need Alen >= "ID", given "ID"\n", need, Alen)) ;
01697   return (FALSE) ;
01698     }
01699
01700     /* since integer overflow has been check, the following cannot overflow: */
01701     Alen -= Col_size + Row_size + (3*n_col + 1) + 5*(n_col+1) + n_row ;
01702
01703     /* Size of A is now Alen >= MAX (2*nnz, 4*n_col) + n_col.  The ordering
01704      * requires Alen >= 2*nnz + n_col, and the postorder requires
01705      * Alen >= 5*n_col. */
01706
01707     ap = Alen ;
01708
01709     /* Front array workspace: 5*(n_col+1) + n_row */
01710     if (!Front_npivcol || !Front_nrows || !Front_ncols || !Front_parent ||
01711         !Front_cols || !Front_cols || !InFront)
01712     {
01713   Front_npivcol = &A [ap] ; ap += (n_col + 1) ;
01714   Front_nrows = &A [ap] ;   ap += (n_col + 1) ;
01715   Front_ncols = &A [ap] ;   ap += (n_col + 1) ;
01716   Front_parent = &A [ap] ;  ap += (n_col + 1) ;
01717   Front_cols = &A [ap] ;    ap += (n_col + 1) ;
01718   InFront = &A [ap] ;   ap += (n_row) ;
01719     }
01720     else
01721     {
01722   /* Fronts are present. Leave the additional space as elbow room. */
01723       ap += 5*(n_col+1) + n_row ;
01724   ap = Alen ;
01725     }
01726
01727     /* Workspace for cset management: 3*n_col+1 */
01728     /* cset_start is of size n_col + 1 */
01729     cset_start = &A [ap] ;
01730     ap += n_col + 1 ;
01731
01732     /* dead_col is of size n_col */
01733     dead_cols = &A [ap] ;
01734     ap += n_col ;
01735
01736     /* cset is of size n_col */
01737     cset = &A [ap] ;
01738     ap += n_col ;
01739
01740     /* Col is of size Col_size.  The space is shared by temp_cstart and csize */
01741     Col = (CColamd_Col *) &A [ap] ;
01742     temp_cstart = (Int *) Col ;   /* [ temp_cstart is of size n_col+1 */
01743     csize = temp_cstart + (n_col+1) ; /* csize is of size n_col+1 */
01744     ap += Col_size ;
01745     ASSERT (Col_size >= 2*n_col+1) ;
01746
01747     /* Row is of size Row_size */
01748     Row = (CColamd_Row *) &A [ap] ;
01749     ap += Row_size ;
01750
01751     /* Initialize csize & dead_cols to zero */
01752     for (i = 0 ; i < n_col ; i++)
01753     {
01754       csize [i] = 0 ;
01755   dead_cols [i] = 0 ;
01756     }
01757
01758     /* === Construct the constraint set ===================================== */
01759
01760     if (n_col == 0)
01761     {
01762   n_cset = 0 ;
01763     }
01764     else if (cmember == (Int *) NULL)
01765     {
01766   /* no constraint set; all columns belong to set zero */
01767   n_cset = 1 ;
01768   csize [0] = n_col ;
01769   DEBUG1 (("no cmember present\n")) ;
01770     }
01771     else
01772     {
01773   n_cset = 0 ;
01774   for (i = 0 ; i < n_col ; i++)
01775   {
01776       if (cmember [i] < 0 || cmember [i] > n_col)
01777       {
01778     stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_invalid_cmember ;
01779     DEBUG1 (("ccolamd: malformed cmember \n")) ;
01780     return (FALSE) ;
01781       }
01782       n_cset = MAX (n_cset, cmember [i]) ;
01783       csize [cmember [i]]++ ;
01784   }
01785   /* cset is zero based */
01786   n_cset++ ;
01787     }
01788
01789     ASSERT ((n_cset >= 0) && (n_cset <= n_col)) ;
01790
01791     cset_start [0] = temp_cstart [0] = 0 ;
01792     for (i = 1 ; i <= n_cset ; i++)
01793     {
01794   cset_start [i] = cset_start [i-1] + csize [i-1] ;
01795   DEBUG4 ((" cset_start ["ID"] = "ID" \n", i , cset_start [i])) ;
01796   temp_cstart [i] = cset_start [i] ;
01797     }
01798
01799     /* do in reverse order to encourage natural tie-breaking */
01800     if (cmember == (Int *) NULL)
01801     {
01802   for (i = n_col-1 ; i >= 0 ; i--)
01803   {
01804       cset [temp_cstart [0]++] = i ;
01805   }
01806     }
01807     else
01808     {
01809   for (i = n_col-1 ; i >= 0 ; i--)
01810   {
01811       cset [temp_cstart [cmember [i]]++] = i ;
01812   }
01813     }
01814
01815     /* ] temp_cstart and csize are no longer used */
01816
01817     /* === Construct the row and column data structures ===================== */
01818
01819     if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats))
01820     {
01821   /* input matrix is invalid */
01822   DEBUG1 (("ccolamd: Matrix invalid\n")) ;
01823   return (FALSE) ;
01824     }
01825
01826     /* === Initialize front info ============================================ */
01827
01828     for (col = 0 ; col < n_col ; col++)
01829     {
01830       Front_npivcol [col] = 0 ;
01831       Front_nrows [col] = 0 ;
01832       Front_ncols [col] = 0 ;
01833       Front_parent [col] = EMPTY ;
01834       Front_cols [col] = EMPTY ;
01835     }
01836
01837     /* === Initialize scores, kill dense rows/columns ======================= */
01838
01839     init_scoring (n_row, n_col, Row, Col, A, p, knobs,
01840   &n_row2, &n_col2, &max_deg, cmember, n_cset, cset_start, dead_cols,
01841   &ndense_row, &nempty_row, &nnewlyempty_row,
01842   &ndense_col, &nempty_col, &nnewlyempty_col) ;
01843
01844     ASSERT (n_row2 == n_row - nempty_row - nnewlyempty_row - ndense_row) ;
01845     ASSERT (n_col2 == n_col - nempty_col - nnewlyempty_col - ndense_col) ;
01846     DEBUG1 (("# dense rows "ID" cols "ID"\n", ndense_row, ndense_col)) ;
01847
01848     /* === Order the supercolumns =========================================== */
01849
01850     ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p,
01851 #ifndef NDEBUG
01852   n_col2,
01853 #endif
01854   max_deg, 2*nnz, cset, cset_start,
01855 #ifndef NDEBUG
01856   n_cset,
01857 #endif
01858   cmember, Front_npivcol, Front_nrows, Front_ncols, Front_parent,
01859   Front_cols, &nfr, aggressive, InFront, order_for_lu) ;
01860
01861     ASSERT (Alen >= 5*n_col) ;
01862
01863     /* === Postorder ======================================================== */
01864
01865     /* A is no longer needed, so use A [0..5*nfr-1] as workspace [ [ */
01866     /* This step requires Alen >= 5*n_col */
01867     Front_child   = A ;
01868     Front_sibling = Front_child + nfr ;
01869     Front_stack   = Front_sibling + nfr ;
01870     Front_order   = Front_stack + nfr ;
01871     Front_size    = Front_order + nfr ;
01872
01873     CCOLAMD_fsize (nfr, Front_size, Front_nrows, Front_ncols,
01874             Front_parent, Front_npivcol) ;
01875
01876     CCOLAMD_postorder (nfr, Front_parent, Front_npivcol, Front_size,
01877         Front_order, Front_child, Front_sibling, Front_stack, Front_cols,
01878   cmember) ;
01879
01880     /* Front_size, Front_stack, Front_child, Front_sibling no longer needed ] */
01881
01882     /* use A [0..nfr-1] as workspace */
01883     CCOLAMD_apply_order (Front_npivcol, Front_order, A, nfr, nfr) ;
01884     CCOLAMD_apply_order (Front_nrows,   Front_order, A, nfr, nfr) ;
01885     CCOLAMD_apply_order (Front_ncols,   Front_order, A, nfr, nfr) ;
01886     CCOLAMD_apply_order (Front_parent,  Front_order, A, nfr, nfr) ;
01887     CCOLAMD_apply_order (Front_cols,    Front_order, A, nfr, nfr) ;
01888
01889     /* fix the parent to refer to the new numbering */
01890     for (i = 0 ; i < nfr ; i++)
01891     {
01892         parent = Front_parent [i] ;
01893         if (parent != EMPTY)
01894         {
01895             Front_parent [i] = Front_order [parent] ;
01896         }
01897     }
01898
01899     /* fix InFront to refer to the new numbering */
01900     for (row = 0 ; row < n_row ; row++)
01901     {
01902         i = InFront [row] ;
01903         ASSERT (i >= EMPTY && i < nfr) ;
01904         if (i != EMPTY)
01905         {
01906             InFront [row] = Front_order [i] ;
01907         }
01908     }
01909
01910     /* Front_order longer needed ] */
01911
01912     /* === Order the columns in the fronts ================================== */
01913
01914     /* use A [0..n_col-1] as inverse permutation */
01915     for (i = 0 ; i < n_col ; i++)
01916     {
01917         A [i] = EMPTY ;
01918     }
01919
01920     k = 0 ;
01921     set1 = 0 ;
01922     for (i = 0 ; i < nfr ; i++)
01923     {
01924         ASSERT (Front_npivcol [i] > 0) ;
01925
01926   set2 = CMEMBER (Front_cols [i]) ;
01927         while (set1 < set2)
01928         {
01929             k += dead_cols [set1] ;
01930             DEBUG3 (("Skip null/dense columns of set "ID"\n",set1)) ;
01931             set1++ ;
01932         }
01933         set1 = set2 ;
01934
01935         for (col = Front_cols [i] ; col != EMPTY ; col = Col [col].nextcol)
01936         {
01937             ASSERT (col >= 0 && col < n_col) ;
01938             DEBUG1 (("ccolamd output ordering: k "ID" col "ID"\n", k, col)) ;
01939             p [k] = col ;
01940             ASSERT (A [col] == EMPTY) ;
01941
01942       cs = CMEMBER (col) ;
01943             ASSERT (k >= cset_start [cs] && k < cset_start [cs+1]) ;
01944
01945             A [col] = k ;
01946             k++ ;
01947         }
01948     }
01949
01950     /* === Order the "dense" and null columns =============================== */
01951
01952     if (n_col2 < n_col)
01953     {
01954         for (col = 0 ; col < n_col ; col++)
01955         {
01956             if (A [col] == EMPTY)
01957             {
01958                 k = Col [col].shared2.order ;
01959     cs = CMEMBER (col) ;
01960 #ifndef NDEBUG
01961                 dead_cols [cs]-- ;
01962 #endif
01963                 ASSERT (k >= cset_start [cs] && k < cset_start [cs+1]) ;
01964                 DEBUG1 (("ccolamd output ordering: k "ID" col "ID
01965                     " (dense or null col)\n", k, col)) ;
01966                 p [k] = col ;
01967                 A [col] = k ;
01968             }
01969         }
01970     }
01971
01972 #ifndef NDEBUG
01973     for (i = 0 ; i < n_cset ; i++)
01974     {
01975       ASSERT (dead_cols [i] == 0) ;
01976     }
01977 #endif
01978
01979     /* === Return statistics in stats ======================================= */
01980
01981     stats [CCOLAMD_DENSE_ROW] = ndense_row ;
01982     stats [CCOLAMD_DENSE_COL] = nempty_row ;
01983     stats [CCOLAMD_NEWLY_EMPTY_ROW] = nnewlyempty_row ;
01984     stats [CCOLAMD_DENSE_COL] = ndense_col ;
01985     stats [CCOLAMD_EMPTY_COL] = nempty_col ;
01986     stats [CCOLAMD_NEWLY_EMPTY_COL] = nnewlyempty_col ;
01987     ASSERT (ndense_col + nempty_col + nnewlyempty_col == n_col - n_col2) ;
01988     if (p_nfr)
01989     {
01990       *p_nfr = nfr ;
01991     }
01992     stats [CCOLAMD_DEFRAG_COUNT] = ngarbage ;
01993     DEBUG1 (("ccolamd: done.\n")) ;
01994     return (TRUE) ;
01995 }
01996
01997
01998 /* ========================================================================== */
01999 /* === colamd_report ======================================================== */
02000 /* ========================================================================== */
02001
02002 PUBLIC void CCOLAMD_report
02003 (
02004     Int stats [CCOLAMD_STATS]
02005 )
02006 {
02007     print_report ("ccolamd", stats) ;
02008 }
02009
02010
02011 /* ========================================================================== */
02012 /* === symamd_report ======================================================== */
02013 /* ========================================================================== */
02014
02015 PUBLIC void CSYMAMD_report
02016 (
02017     Int stats [CCOLAMD_STATS]
02018 )
02019 {
02020     print_report ("csymamd", stats) ;
02021 }
02022
02023
02024 /* ========================================================================== */
02025 /* === NON-USER-CALLABLE ROUTINES: ========================================== */
02026 /* ========================================================================== */
02027
02028 /* There are no user-callable routines beyond this point in the file */
02029
02030
02031 /* ========================================================================== */
02032 /* === init_rows_cols ======================================================= */
02033 /* ========================================================================== */
02034
02035 /*
02036     Takes the column form of the matrix in A and creates the row form of the
02037     matrix.  Also, row and column attributes are stored in the Col and Row
02038     structs.  If the columns are un-sorted or contain duplicate row indices,
02039     this routine will also sort and remove duplicate row indices from the
02040     column form of the matrix.  Returns FALSE if the matrix is invalid,
02041     TRUE otherwise.  Not user-callable.
02042 */
02043
02044 PRIVATE Int init_rows_cols  /* returns TRUE if OK, or FALSE otherwise */
02045 (
02046     /* === Parameters ======================================================= */
02047
02048     Int n_row,      /* number of rows of A */
02049     Int n_col,      /* number of columns of A */
02050     CColamd_Row Row [ ],    /* of size n_row+1 */
02051     CColamd_Col Col [ ],    /* of size n_col+1 */
02052     Int A [ ],      /* row indices of A, of size Alen */
02053     Int p [ ],      /* pointers to columns in A, of size n_col+1 */
02054     Int stats [CCOLAMD_STATS] /* colamd statistics */
02055 )
02056 {
02057     /* === Local variables ================================================== */
02058
02059     Int col ;     /* a column index */
02060     Int row ;     /* a row index */
02061     Int *cp ;     /* a column pointer */
02062     Int *cp_end ;   /* a pointer to the end of a column */
02063     Int *rp ;     /* a row pointer */
02064     Int *rp_end ;   /* a pointer to the end of a row */
02065     Int last_row ;    /* previous row */
02066
02067     /* === Initialize columns, and check column pointers ==================== */
02068
02069     for (col = 0 ; col < n_col ; col++)
02070     {
02071   Col [col].start = p [col] ;
02072   Col [col].length = p [col+1] - p [col] ;
02073
02074   if (Col [col].length < 0)
02075   {
02076       /* column pointers must be non-decreasing */
02077       stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_col_length_negative ;
02078       stats [CCOLAMD_INFO1] = col ;
02079       stats [CCOLAMD_INFO2] = Col [col].length ;
02080       DEBUG1 (("ccolamd: col "ID" length "ID" < 0\n",
02081       col, Col [col].length)) ;
02082       return (FALSE) ;
02083   }
02084
02085   Col [col].shared1.thickness = 1 ;
02086   Col [col].shared2.score = 0 ;
02087   Col [col].shared3.prev = EMPTY ;
02088   Col [col].shared4.degree_next = EMPTY ;
02089         Col [col].nextcol = EMPTY ;
02090         Col [col].lastcol = col ;
02091     }
02092
02093     /* p [0..n_col] no longer needed, used as "head" in subsequent routines */
02094
02095     /* === Scan columns, compute row degrees, and check row indices ========= */
02096
02097     stats [CCOLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/
02098
02099     for (row = 0 ; row < n_row ; row++)
02100     {
02101   Row [row].length = 0 ;
02102   Row [row].shared2.mark = -1 ;
02103         Row [row].thickness = 1 ;
02104         Row [row].front = EMPTY ;
02105     }
02106
02107     for (col = 0 ; col < n_col ; col++)
02108     {
02109   DEBUG1 (("\nCcolamd input column "ID":\n", col)) ;
02110   last_row = -1 ;
02111
02112   cp = &A [p [col]] ;
02113   cp_end = &A [p [col+1]] ;
02114
02115   while (cp < cp_end)
02116   {
02117       row = *cp++ ;
02118       DEBUG1 (("row: "ID"\n", row)) ;
02119
02120       /* make sure row indices within range */
02121       if (row < 0 || row >= n_row)
02122       {
02123     stats [CCOLAMD_STATUS] = CCOLAMD_ERROR_row_index_out_of_bounds ;
02124     stats [CCOLAMD_INFO1] = col ;
02125     stats [CCOLAMD_INFO2] = row ;
02126     stats [CCOLAMD_INFO3] = n_row ;
02127     DEBUG1 (("row "ID" col "ID" out of bounds\n", row, col)) ;
02128     return (FALSE) ;
02129       }
02130
02131       if (row <= last_row || Row [row].shared2.mark == col)
02132       {
02133     /* row index are unsorted or repeated (or both), thus col */
02134     /* is jumbled.  This is a notice, not an error condition. */
02135     stats [CCOLAMD_STATUS] = CCOLAMD_OK_BUT_JUMBLED ;
02136     stats [CCOLAMD_INFO1] = col ;
02137     stats [CCOLAMD_INFO2] = row ;
02138     (stats [CCOLAMD_INFO3]) ++ ;
02139     DEBUG1 (("row "ID" col "ID" unsorted/duplicate\n", row, col)) ;
02140       }
02141
02142       if (Row [row].shared2.mark != col)
02143       {
02144     Row [row].length++ ;
02145       }
02146       else
02147       {
02148     /* this is a repeated entry in the column, */
02149     /* it will be removed */
02150     Col [col].length-- ;
02151       }
02152
02153       /* mark the row as having been seen in this column */
02154       Row [row].shared2.mark = col ;
02155
02156       last_row = row ;
02157   }
02158     }
02159
02160     /* === Compute row pointers ============================================= */
02161
02162     /* row form of the matrix starts directly after the column */
02163     /* form of matrix in A */
02164     Row [0].start = p [n_col] ;
02165     Row [0].shared1.p = Row [0].start ;
02166     Row [0].shared2.mark = -1 ;
02167     for (row = 1 ; row < n_row ; row++)
02168     {
02169   Row [row].start = Row [row-1].start + Row [row-1].length ;
02170   Row [row].shared1.p = Row [row].start ;
02171   Row [row].shared2.mark = -1 ;
02172     }
02173
02174     /* === Create row form ================================================== */
02175
02176     if (stats [CCOLAMD_STATUS] == CCOLAMD_OK_BUT_JUMBLED)
02177     {
02178   /* if cols jumbled, watch for repeated row indices */
02179   for (col = 0 ; col < n_col ; col++)
02180   {
02181       cp = &A [p [col]] ;
02182       cp_end = &A [p [col+1]] ;
02183       while (cp < cp_end)
02184       {
02185     row = *cp++ ;
02186     if (Row [row].shared2.mark != col)
02187     {
02188         A [(Row [row].shared1.p)++] = col ;
02189         Row [row].shared2.mark = col ;
02190     }
02191       }
02192   }
02193     }
02194     else
02195     {
02196   /* if cols not jumbled, we don't need the mark (this is faster) */
02197   for (col = 0 ; col < n_col ; col++)
02198   {
02199       cp = &A [p [col]] ;
02200       cp_end = &A [p [col+1]] ;
02201       while (cp < cp_end)
02202       {
02203     A [(Row [*cp++].shared1.p)++] = col ;
02204       }
02205   }
02206     }
02207
02208     /* === Clear the row marks and set row degrees ========================== */
02209
02210     for (row = 0 ; row < n_row ; row++)
02211     {
02212   Row [row].shared2.mark = 0 ;
02213   Row [row].shared1.degree = Row [row].length ;
02214     }
02215
02216     /* === See if we need to re-create columns ============================== */
02217
02218     if (stats [CCOLAMD_STATUS] == CCOLAMD_OK_BUT_JUMBLED)
02219     {
02220       DEBUG1 (("ccolamd: reconstructing column form, matrix jumbled\n")) ;
02221
02222 #ifndef NDEBUG
02223   /* make sure column lengths are correct */
02224   for (col = 0 ; col < n_col ; col++)
02225   {
02226       p [col] = Col [col].length ;
02227   }
02228   for (row = 0 ; row < n_row ; row++)
02229   {
02230       rp = &A [Row [row].start] ;
02231       rp_end = rp + Row [row].length ;
02232       while (rp < rp_end)
02233       {
02234     p [*rp++]-- ;
02235       }
02236   }
02237   for (col = 0 ; col < n_col ; col++)
02238   {
02239       ASSERT (p [col] == 0) ;
02240   }
02241   /* now p is all zero (different than when debugging is turned off) */
02242 #endif
02243
02244   /* === Compute col pointers ========================================= */
02245
02246   /* col form of the matrix starts at A [0]. */
02247   /* Note, we may have a gap between the col form and the row */
02248   /* form if there were duplicate entries, if so, it will be */
02249   /* removed upon the first garbage collection */
02250   Col [0].start = 0 ;
02251   p [0] = Col [0].start ;
02252   for (col = 1 ; col < n_col ; col++)
02253   {
02254       /* note that the lengths here are for pruned columns, i.e. */
02255       /* no duplicate row indices will exist for these columns */
02256       Col [col].start = Col [col-1].start + Col [col-1].length ;
02257       p [col] = Col [col].start ;
02258   }
02259
02260   /* === Re-create col form =========================================== */
02261
02262   for (row = 0 ; row < n_row ; row++)
02263   {
02264       rp = &A [Row [row].start] ;
02265       rp_end = rp + Row [row].length ;
02266       while (rp < rp_end)
02267       {
02268     A [(p [*rp++])++] = row ;
02269       }
02270   }
02271     }
02272
02273     /* === Done.  Matrix is not (or no longer) jumbled ====================== */
02274
02275
02276     return (TRUE) ;
02277 }
02278
02279
02280 /* ========================================================================== */
02281 /* === init_scoring ========================================================= */
02282 /* ========================================================================== */
02283
02284 /*
02285     Kills dense or empty columns and rows, calculates an initial score for
02286     each column, and places all columns in the degree lists.  Not user-callable.
02287 */
02288
02289 PRIVATE void init_scoring
02290 (
02291     /* === Parameters ======================================================= */
02292
02293     Int n_row,      /* number of rows of A */
02294     Int n_col,      /* number of columns of A */
02295     CColamd_Row Row [ ],  /* of size n_row+1 */
02296     CColamd_Col Col [ ],  /* of size n_col+1 */
02297     Int A [ ],      /* column form and row form of A */
02298     Int head [ ],   /* of size n_col+1 */
02299     double knobs [CCOLAMD_KNOBS],/* parameters */
02300     Int *p_n_row2,    /* number of non-dense, non-empty rows */
02301     Int *p_n_col2,    /* number of non-dense, non-empty columns */
02302     Int *p_max_deg,   /* maximum row degree */
02303     Int cmember [ ],
02304     Int n_cset,
02305     Int cset_start [ ],
02306     Int dead_cols [ ],
02307     Int *p_ndense_row,    /* number of dense rows */
02308     Int *p_nempty_row,    /* number of original empty rows */
02309     Int *p_nnewlyempty_row, /* number of newly empty rows */
02310     Int *p_ndense_col,    /* number of dense cols (excl "empty" cols) */
02311     Int *p_nempty_col,    /* number of original empty cols */
02312     Int *p_nnewlyempty_col  /* number of newly empty cols */
02313 )
02314 {
02315 /* === Local variables ================================================== */
02316
02317     Int c ;     /* a column index */
02318     Int r, row ;    /* a row index */
02319     Int *cp ;     /* a column pointer */
02320     Int deg ;     /* degree of a row or column */
02321     Int *cp_end ;   /* a pointer to the end of a column */
02322     Int *new_cp ;   /* new column pointer */
02323     Int col_length ;    /* length of pruned column */
02324     Int score ;     /* current column score */
02325     Int n_col2 ;    /* number of non-dense, non-empty columns */
02326     Int n_row2 ;    /* number of non-dense, non-empty rows */
02327     Int dense_row_count ; /* remove rows with more entries than this */
02328     Int dense_col_count ; /* remove cols with more entries than this */
02329     Int max_deg ;   /* maximum row degree */
02330     Int s ;     /* a cset index */
02331     Int ndense_row ;    /* number of dense rows */
02332     Int nempty_row ;    /* number of empty rows */
02333     Int nnewlyempty_row ; /* number of newly empty rows */
02334     Int ndense_col ;    /* number of dense cols (excl "empty" cols) */
02335     Int nempty_col ;    /* number of original empty cols */
02336     Int nnewlyempty_col ; /* number of newly empty cols */
02337     Int ne ;
02338
02339 #ifndef NDEBUG
02340     Int debug_count ;   /* debug only. */
02341 #endif
02342
02343     /* === Extract knobs ==================================================== */
02344
02345     /* Note: if knobs contains a NaN, this is undefined: */
02346     if (knobs [CCOLAMD_DENSE_ROW] < 0)
02347     {
02348   /* only remove completely dense rows */
02349   dense_row_count = n_col-1 ;
02350     }
02351     else
02352     {
02353   dense_row_count = DENSE_DEGREE (knobs [CCOLAMD_DENSE_ROW], n_col) ;
02354     }
02355     if (knobs [CCOLAMD_DENSE_COL] < 0)
02356     {
02357   /* only remove completely dense columns */
02358   dense_col_count = n_row-1 ;
02359     }
02360     else
02361     {
02362   dense_col_count =
02363       DENSE_DEGREE (knobs [CCOLAMD_DENSE_COL], MIN (n_row, n_col)) ;
02364     }
02365
02366     DEBUG1 (("densecount: "ID" "ID"\n", dense_row_count, dense_col_count)) ;
02367     max_deg = 0 ;
02368
02369     n_col2 = n_col ;
02370     n_row2 = n_row ;
02371
02372     /* Set the head array for bookkeeping of dense and empty columns. */
02373     /* This will be used as hash buckets later. */
02374     for (s = 0 ; s < n_cset ; s++)
02375     {
02376   head [s] = cset_start [s+1] ;
02377     }
02378
02379     ndense_col = 0 ;
02380     nempty_col = 0 ;
02381     nnewlyempty_col = 0 ;
02382     ndense_row = 0 ;
02383     nempty_row = 0 ;
02384     nnewlyempty_row = 0 ;
02385
02386     /* === Kill empty columns =============================================== */
02387
02388     /* Put the empty columns at the end in their natural order, so that LU */
02389     /* factorization can proceed as far as possible. */
02390     for (c = n_col-1 ; c >= 0 ; c--)
02391     {
02392   deg = Col [c].length ;
02393   if (deg == 0)
02394   {
02395       /* this is a empty column, kill and order it last of its cset */
02396       Col [c].shared2.order = --head [CMEMBER (c)] ;
02397       --n_col2 ;
02398       dead_cols [CMEMBER (c)] ++ ;
02399       nempty_col++ ;
02400       KILL_PRINCIPAL_COL (c) ;
02401   }
02402     }
02403     DEBUG1 (("ccolamd: null columns killed: "ID"\n", n_col - n_col2)) ;
02404
02405     /* === Kill dense columns =============================================== */
02406
02407     /* Put the dense columns at the end, in their natural order */
02408     for (c = n_col-1 ; c >= 0 ; c--)
02409     {
02410   /* skip any dead columns */
02411   if (COL_IS_DEAD (c))
02412   {
02413       continue ;
02414   }
02415   deg = Col [c].length ;
02416   if (deg > dense_col_count)
02417   {
02418       /* this is a dense column, kill and order it last of its cset */
02419       Col [c].shared2.order = --head [CMEMBER (c)] ;
02420       --n_col2 ;
02421       dead_cols [CMEMBER (c)] ++ ;
02422       ndense_col++ ;
02423       /* decrement the row degrees */
02424       cp = &A [Col [c].start] ;
02425       cp_end = cp + Col [c].length ;
02426       while (cp < cp_end)
02427       {
02428     Row [*cp++].shared1.degree-- ;
02429       }
02430       KILL_PRINCIPAL_COL (c) ;
02431   }
02432     }
02433     DEBUG1 (("Dense and null columns killed: "ID"\n", n_col - n_col2)) ;
02434
02435     /* === Kill dense and empty rows ======================================== */
02436
02437     /* Note that there can now be empty rows, since dense columns have
02438      * been deleted.  These are "newly" empty rows. */
02439
02440     ne = 0 ;
02441     for (r = 0 ; r < n_row ; r++)
02442     {
02443   deg = Row [r].shared1.degree ;
02444   ASSERT (deg >= 0 && deg <= n_col) ;
02445         if (deg > dense_row_count)
02446         {
02447             /* There is at least one dense row.  Continue ordering, but */
02448             /* symbolic factorization will be redone after ccolamd is done.*/
02449             ndense_row++ ;
02450         }
02451         if (deg == 0)
02452         {
02453             /* this is a newly empty row, or original empty row */
02454             ne++ ;
02455         }
02456   if (deg > dense_row_count || deg == 0)
02457   {
02458       /* kill a dense or empty row */
02459       KILL_ROW (r) ;
02460       Row [r].thickness = 0 ;
02461       --n_row2 ;
02462   }
02463   else
02464   {
02465       /* keep track of max degree of remaining rows */
02466       max_deg = MAX (max_deg, deg) ;
02467   }
02468     }
02469     nnewlyempty_row = ne - nempty_row ;
02470     DEBUG1 (("ccolamd: Dense and null rows killed: "ID"\n", n_row - n_row2)) ;
02471
02472     /* === Compute initial column scores ==================================== */
02473
02474     /* At this point the row degrees are accurate.  They reflect the number */
02475     /* of "live" (non-dense) columns in each row.  No empty rows exist. */
02476     /* Some "live" columns may contain only dead rows, however.  These are */
02477     /* pruned in the code below. */
02478
02479     /* now find the initial COLMMD score for each column */
02480     for (c = n_col-1 ; c >= 0 ; c--)
02481     {
02482   /* skip dead column */
02483   if (COL_IS_DEAD (c))
02484   {
02485       continue ;
02486   }
02487   score = 0 ;
02488   cp = &A [Col [c].start] ;
02489   new_cp = cp ;
02490   cp_end = cp + Col [c].length ;
02491   while (cp < cp_end)
02492   {
02493       /* get a row */
02494       row = *cp++ ;
02495       /* skip if dead */
02496       if (ROW_IS_DEAD (row))
02497       {
02498     continue ;
02499       }
02500       /* compact the column */
02501       *new_cp++ = row ;
02502       /* add row's external degree */
02503       score += Row [row].shared1.degree - 1 ;
02504       /* guard against integer overflow */
02505       score = MIN (score, n_col) ;
02506   }
02507   /* determine pruned column length */
02508   col_length = (Int) (new_cp - &A [Col [c].start]) ;
02509   if (col_length == 0)
02510   {
02511       /* a newly-made null column (all rows in this col are "dense" */
02512       /* and have already been killed) */
02513       DEBUG1 (("Newly null killed: "ID"\n", c)) ;
02514       Col [c].shared2.order = -- head [CMEMBER (c)] ;
02515       --n_col2 ;
02516       dead_cols [CMEMBER (c)] ++ ;
02517             nnewlyempty_col++ ;
02518       KILL_PRINCIPAL_COL (c) ;
02519   }
02520   else
02521   {
02522       /* set column length and set score */
02523       ASSERT (score >= 0) ;
02524       ASSERT (score <= n_col) ;
02525       Col [c].length = col_length ;
02526       Col [c].shared2.score = score ;
02527   }
02528     }
02529     DEBUG1 (("ccolamd: Dense, null, and newly-null columns killed: "ID"\n",
02530       n_col-n_col2)) ;
02531
02532     /* At this point, all empty rows and columns are dead.  All live columns */
02533     /* are "clean" (containing no dead rows) and simplicial (no supercolumns */
02534     /* yet).  Rows may contain dead columns, but all live rows contain at */
02535     /* least one live column. */
02536
02537 #ifndef NDEBUG
02538     debug_count = 0 ;
02539 #endif
02540
02541     /* clear the hash buckets */
02542     for (c = 0 ; c <= n_col ; c++)
02543     {
02544   head [c] = EMPTY ;
02545     }
02546
02547 #ifndef NDEBUG
02548     debug_structures (n_row, n_col, Row, Col, A, cmember, cset_start) ;
02549 #endif
02550
02551     /* === Return number of remaining columns, and max row degree =========== */
02552
02553     *p_n_col2 = n_col2 ;
02554     *p_n_row2 = n_row2 ;
02555     *p_max_deg = max_deg ;
02556     *p_ndense_row = ndense_row ;
02557     *p_nempty_row = nempty_row ;        /* original empty rows */
02558     *p_nnewlyempty_row = nnewlyempty_row ;
02559     *p_ndense_col = ndense_col ;
02560     *p_nempty_col = nempty_col ;        /* original empty cols */
02561     *p_nnewlyempty_col = nnewlyempty_col ;
02562 }
02563
02564
02565 /* ========================================================================== */
02566 /* === find_ordering ======================================================== */
02567 /* ========================================================================== */
02568
02569 /*
02570  *   Order the principal columns of the supercolumn form of the matrix
02571  *  (no supercolumns on input).  Uses a minimum approximate column minimum
02572  *  degree ordering method.  Not user-callable.
02573  */
02574
02575 PRIVATE Int find_ordering /* return the number of garbage collections */
02576 (
02577     /* === Parameters ======================================================= */
02578
02579     Int n_row,      /* number of rows of A */
02580     Int n_col,      /* number of columns of A */
02581     Int Alen,     /* size of A, 2*nnz + n_col or larger */
02582     CColamd_Row Row [ ],  /* of size n_row+1 */
02583     CColamd_Col Col [ ],  /* of size n_col+1 */
02584     Int A [ ],      /* column form and row form of A */
02585     Int head [ ],   /* of size n_col+1 */
02586 #ifndef NDEBUG
02587     Int n_col2,     /* Remaining columns to order */
02588 #endif
02589     Int max_deg,    /* Maximum row degree */
02590     Int pfree,      /* index of first free slot (2*nnz on entry) */
02591     Int cset [ ],   /* constraint set of A */
02592     Int cset_start [ ],   /* pointer to the start of every cset */
02593 #ifndef NDEBUG
02594     Int n_cset,     /* number of csets */
02595 #endif
02596     Int cmember [ ],    /* col -> cset mapping */
02597     Int Front_npivcol [ ],
02598     Int Front_nrows [ ],
02599     Int Front_ncols [ ],
02600     Int Front_parent [ ],
02601     Int Front_cols [ ],
02602     Int *p_nfr,                /* number of fronts */
02603     Int aggressive,
02604     Int InFront [ ],
02605     Int order_for_lu
02606 )
02607 {
02608     /* === Local variables ================================================== */
02609
02610     Int k ;     /* current pivot ordering step */
02611     Int pivot_col ;   /* current pivot column */
02612     Int *cp ;     /* a column pointer */
02613     Int *rp ;     /* a row pointer */
02614     Int pivot_row ;   /* current pivot row */
02615     Int *new_cp ;   /* modified column pointer */
02616     Int *new_rp ;   /* modified row pointer */
02617     Int pivot_row_start ; /* pointer to start of pivot row */
02618     Int pivot_row_degree ;  /* number of columns in pivot row */
02619     Int pivot_row_length ;  /* number of supercolumns in pivot row */
02620     Int pivot_col_score ; /* score of pivot column */
02621     Int needed_memory ;   /* free space needed for pivot row */
02622     Int *cp_end ;   /* pointer to the end of a column */
02623     Int *rp_end ;   /* pointer to the end of a row */
02624     Int row ;     /* a row index */
02625     Int col ;     /* a column index */
02626     Int max_score ;   /* maximum possible score */
02627     Int cur_score ;   /* score of current column */
02628     unsigned Int hash ;   /* hash value for supernode detection */
02629     Int head_column ;   /* head of hash bucket */
02630     Int first_col ;   /* first column in hash bucket */
02631     Int tag_mark ;    /* marker value for mark array */
02632     Int row_mark ;    /* Row [row].shared2.mark */
02633     Int set_difference ;  /* set difference size of row with pivot row */
02634     Int min_score ;   /* smallest column score */
02635     Int col_thickness ;   /* "thickness" (no. of columns in a supercol) */
02636     Int max_mark ;    /* maximum value of tag_mark */
02637     Int pivot_col_thickness ; /* number of columns represented by pivot col */
02638     Int prev_col ;    /* Used by Dlist operations. */
02639     Int next_col ;    /* Used by Dlist operations. */
02640     Int ngarbage ;    /* number of garbage collections performed */
02641     Int current_set ;   /* consraint set that is being ordered */
02642     Int score ;     /* score of a column */
02643     Int colstart ;    /* pointer to first column in current cset */
02644     Int colend ;    /* pointer to last column in current cset */
02645     Int deadcol ;   /* number of dense & null columns in a cset */
02646
02647 #ifndef NDEBUG
02648     Int debug_d ;   /* debug loop counter */
02649     Int debug_step = 0 ;  /* debug loop counter */
02650     Int cols_thickness = 0 ;  /* the thickness of the columns in current */
02651             /* cset degreelist and in pivot row pattern. */
02652 #endif
02653
02654     Int pivot_row_thickness ;   /* number of rows represented by pivot row */
02655     Int nfr = 0 ;               /* number of fronts */
02656     Int child ;
02657
02658     /* === Initialization and clear mark ==================================== */
02659
02660     max_mark = Int_MAX - n_col ;  /* Int_MAX defined in <limits.h> */
02661     tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02662     min_score = 0 ;
02663     ngarbage = 0 ;
02664     current_set = -1 ;
02665     deadcol = 0 ;
02666     DEBUG1 (("ccolamd: Ordering, n_col2="ID"\n", n_col2)) ;
02667
02668     for (row = 0 ; row < n_row ; row++)
02669     {
02670         InFront [row] = EMPTY ;
02671     }
02672
02673     /* === Order the columns ================================================ */
02674
02675     for (k = 0 ; k < n_col ; /* 'k' is incremented below */)
02676     {
02677
02678   /* make sure degree list isn't empty */
02679   ASSERT (min_score >= 0) ;
02680   ASSERT (min_score <= n_col) ;
02681   ASSERT (head [min_score] >= EMPTY) ;
02682
02683 #ifndef NDEBUG
02684   for (debug_d = 0 ; debug_d < min_score ; debug_d++)
02685   {
02686       ASSERT (head [debug_d] == EMPTY) ;
02687   }
02688 #endif
02689
02690   /* Initialize the degree list with columns from next non-empty cset */
02691
02692   while ((k+deadcol) == cset_start [current_set+1])
02693   {
02694       current_set++ ;
02695       DEBUG1 (("\n\n\n============ CSET: "ID"\n", current_set)) ;
02696       k += deadcol ;  /* jump to start of next cset */
02697         deadcol = 0 ; /* reset dead column count */
02698
02699       ASSERT ((current_set == n_cset) == (k == n_col)) ;
02700
02701       /* return if all columns are ordered. */
02702       if (k == n_col)
02703       {
02704     *p_nfr = nfr ;
02705         return (ngarbage) ;
02706       }
02707
02708 #ifndef NDEBUG
02709       for (col = 0 ; col <= n_col ; col++)
02710       {
02711           ASSERT (head [col] == EMPTY) ;
02712       }
02713 #endif
02714
02715       min_score = n_col ;
02716       colstart = cset_start [current_set] ;
02717       colend = cset_start [current_set+1] ;
02718
02719       while (colstart < colend)
02720       {
02721     col = cset [colstart++] ;
02722
02723     if (COL_IS_DEAD(col))
02724     {
02725         DEBUG1 (("Column "ID" is dead\n", col)) ;
02726         /* count dense and null columns */
02727         if (Col [col].shared2.order != EMPTY)
02728         {
02729       deadcol++ ;
02730         }
02731         continue ;
02732     }
02733
02734     /* only add principal columns in current set to degree lists */
02735     ASSERT (CMEMBER (col) == current_set) ;
02736
02737     score = Col [col].shared2.score ;
02738     DEBUG1 (("Column "ID" is alive, score "ID"\n", col, score)) ;
02739
02740     ASSERT (min_score >= 0) ;
02741     ASSERT (min_score <= n_col) ;
02742     ASSERT (score >= 0) ;
02743     ASSERT (score <= n_col) ;
02744     ASSERT (head [score] >= EMPTY) ;
02745
02746     /* now add this column to dList at proper score location */
02747     next_col = head [score] ;
02748     Col [col].shared3.prev = EMPTY ;
02749     Col [col].shared4.degree_next = next_col ;
02750
02751     /* if there already was a column with the same score, set its */
02752     /* previous pointer to this new column */
02753     if (next_col != EMPTY)
02754     {
02755         Col [next_col].shared3.prev = col ;
02756     }
02757     head [score] = col ;
02758
02759     /* see if this score is less than current min */
02760     min_score = MIN (min_score, score) ;
02761       }
02762
02763 #ifndef NDEBUG
02764       DEBUG1 (("degree lists initialized \n")) ;
02765       debug_deg_lists (n_row, n_col, Row, Col, head, min_score,
02766     ((cset_start [current_set+1]-cset_start [current_set])-deadcol),
02767     max_deg) ;
02768 #endif
02769   }
02770
02771 #ifndef NDEBUG
02772   if (debug_step % 100 == 0)
02773   {
02774       DEBUG2 (("\n...   Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
02775   }
02776   else
02777   {
02778       DEBUG3 (("\n------Step k: "ID" out of n_col2: "ID"\n", k, n_col2)) ;
02779   }
02780   debug_step++ ;
02781   DEBUG1 (("start of step k="ID": ", k)) ;
02782   debug_deg_lists (n_row, n_col, Row, Col, head,
02783        min_score, cset_start [current_set+1]-(k+deadcol), max_deg) ;
02784   debug_matrix (n_row, n_col, Row, Col, A) ;
02785 #endif
02786
02787   /* === Select pivot column, and order it ============================ */
02788
02789   while (head [min_score] == EMPTY && min_score < n_col)
02790   {
02791       min_score++ ;
02792   }
02793
02794   pivot_col = head [min_score] ;
02795
02796   ASSERT (pivot_col >= 0 && pivot_col <= n_col) ;
02797   next_col = Col [pivot_col].shared4.degree_next ;
02798   head [min_score] = next_col ;
02799   if (next_col != EMPTY)
02800   {
02801       Col [next_col].shared3.prev = EMPTY ;
02802   }
02803
02804   ASSERT (COL_IS_ALIVE (pivot_col)) ;
02805
02806   /* remember score for defrag check */
02807   pivot_col_score = Col [pivot_col].shared2.score ;
02808
02809   /* the pivot column is the kth column in the pivot order */
02810   Col [pivot_col].shared2.order = k ;
02811
02812   /* increment order count by column thickness */
02813   pivot_col_thickness = Col [pivot_col].shared1.thickness ;
02814   k += pivot_col_thickness ;
02815   ASSERT (pivot_col_thickness > 0) ;
02816   DEBUG3 (("Pivot col: "ID" thick "ID"\n", pivot_col,
02817         pivot_col_thickness)) ;
02818
02819   /* === Garbage_collection, if necessary ============================= */
02820
02821   needed_memory = MIN (pivot_col_score, n_col - k) ;
02822   if (pfree + needed_memory >= Alen)
02823   {
02824       pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ;
02825       ngarbage++ ;
02826       /* after garbage collection we will have enough */
02827       ASSERT (pfree + needed_memory < Alen) ;
02828       /* garbage collection has wiped out Row [ ].shared2.mark array */
02829       tag_mark = clear_mark (0, max_mark, n_row, Row) ;
02830
02831 #ifndef NDEBUG
02832       debug_matrix (n_row, n_col, Row, Col, A) ;
02833 #endif
02834   }
02835
02836   /* === Compute pivot row pattern ==================================== */
02837
02838   /* get starting location for this new merged row */
02839   pivot_row_start = pfree ;
02840
02841   /* initialize new row counts to zero */
02842   pivot_row_degree = 0 ;
02843         pivot_row_thickness = 0 ;
02844
02845   /* tag pivot column as having been visited so it isn't included */
02846   /* in merged pivot row */
02847   Col [pivot_col].shared1.thickness = -pivot_col_thickness ;
02848
02849   /* pivot row is the union of all rows in the pivot column pattern */
02850   cp = &A [Col [pivot_col].start] ;
02851   cp_end = cp + Col [pivot_col].length ;
02852   while (cp < cp_end)
02853   {
02854       /* get a row */
02855       row = *cp++ ;
02856       ASSERT (row >= 0 && row < n_row) ;
02857       DEBUG4 (("Pivcol pattern "ID" "ID"\n", ROW_IS_ALIVE (row), row)) ;
02858       /* skip if row is dead */
02859       if (ROW_IS_ALIVE (row))
02860       {
02861     /* sum the thicknesses of all the rows */
02862     pivot_row_thickness += Row [row].thickness ;
02863
02864     rp = &A [Row [row].start] ;
02865     rp_end = rp + Row [row].length ;
02866     while (rp < rp_end)
02867     {
02868         /* get a column */
02869         col = *rp++ ;
02870         /* add the column, if alive and untagged */
02871         col_thickness = Col [col].shared1.thickness ;
02872         if (col_thickness > 0 && COL_IS_ALIVE (col))
02873         {
02874       /* tag column in pivot row */
02875       Col [col].shared1.thickness = -col_thickness ;
02876       ASSERT (pfree < Alen) ;
02877       /* place column in pivot row */
02878       A [pfree++] = col ;
02879       pivot_row_degree += col_thickness ;
02880       DEBUG4 (("\t\t\tNew live col in pivrow: "ID"\n",col)) ;
02881         }
02882 #ifndef NDEBUG
02883         if (col_thickness < 0 && COL_IS_ALIVE (col))
02884         {
02885       DEBUG4 (("\t\t\tOld live col in pivrow: "ID"\n",col)) ;
02886         }
02887 #endif
02888     }
02889       }
02890   }
02891
02892         /* pivot_row_thickness is the number of rows in frontal matrix */
02893         /* including both pivotal rows and nonpivotal rows */
02894
02895   /* clear tag on pivot column */
02896   Col [pivot_col].shared1.thickness = pivot_col_thickness ;
02897   max_deg = MAX (max_deg, pivot_row_degree) ;
02898
02899 #ifndef NDEBUG
02900   DEBUG3 (("check2\n")) ;
02901   debug_mark (n_row, Row, tag_mark, max_mark) ;
02902 #endif
02903
02904   /* === Kill all rows used to construct pivot row ==================== */
02905
02906   /* also kill pivot row, temporarily */
02907   cp = &A [Col [pivot_col].start] ;
02908   cp_end = cp + Col [pivot_col].length ;
02909   while (cp < cp_end)
02910   {
02911       /* may be killing an already dead row */
02912       row = *cp++ ;
02913       DEBUG3 (("Kill row in pivot col: "ID"\n", row)) ;
02914       ASSERT (row >= 0 && row < n_row) ;
02915             if (ROW_IS_ALIVE (row))
02916             {
02917                 if (Row [row].front != EMPTY)
02918                 {
02919                     /* This row represents a frontal matrix. */
02920                     /* Row [row].front is a child of current front */
02921                     child = Row [row].front ;
02922                     Front_parent [child] = nfr ;
02923                     DEBUG1 (("Front "ID" => front "ID", normal\n", child, nfr));
02924                 }
02925                 else
02926                 {
02927                     /* This is an original row.  Keep track of which front
02928                      * is its parent in the row-merge tree. */
02929                     InFront [row] = nfr ;
02930                     DEBUG1 (("Row "ID" => front "ID", normal\n", row, nfr)) ;
02931                 }
02932             }
02933
02934             KILL_ROW (row) ;
02935             Row [row].thickness = 0 ;
02936   }
02937
02938   /* === Select a row index to use as the new pivot row =============== */
02939
02940   pivot_row_length = pfree - pivot_row_start ;
02941   if (pivot_row_length > 0)
02942   {
02943       /* pick the "pivot" row arbitrarily (first row in col) */
02944       pivot_row = A [Col [pivot_col].start] ;
02945       DEBUG3 (("Pivotal row is "ID"\n", pivot_row)) ;
02946   }
02947   else
02948   {
02949       /* there is no pivot row, since it is of zero length */
02950       pivot_row = EMPTY ;
02951       ASSERT (pivot_row_length == 0) ;
02952   }
02953   ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ;
02954
02955   /* === Approximate degree computation =============================== */
02956
02957   /* Here begins the computation of the approximate degree.  The column */
02958   /* score is the sum of the pivot row "length", plus the size of the */
02959   /* set differences of each row in the column minus the pattern of the */
02960   /* pivot row itself.  The column ("thickness") itself is also */
02961   /* excluded from the column score (we thus use an approximate */
02962   /* external degree). */
02963
02964   /* The time taken by the following code (compute set differences, and */
02965   /* add them up) is proportional to the size of the data structure */
02966   /* being scanned - that is, the sum of the sizes of each column in */
02967   /* the pivot row.  Thus, the amortized time to compute a column score */
02968   /* is proportional to the size of that column (where size, in this */
02969   /* context, is the column "length", or the number of row indices */
02970   /* in that column).  The number of row indices in a column is */
02971   /* monotonically non-decreasing, from the length of the original */
02972   /* column on input to colamd. */
02973
02974   /* === Compute set differences ====================================== */
02975
02976   DEBUG3 (("** Computing set differences phase. **\n")) ;
02977
02978   /* pivot row is currently dead - it will be revived later. */
02979
02980   DEBUG3 (("Pivot row: ")) ;
02981   /* for each column in pivot row */
02982   rp = &A [pivot_row_start] ;
02983   rp_end = rp + pivot_row_length ;
02984   while (rp < rp_end)
02985   {
02986       col = *rp++ ;
02987       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
02988       DEBUG3 (("Col: "ID"\n", col)) ;
02989
02990       /* clear tags used to construct pivot row pattern */
02991       col_thickness = -Col [col].shared1.thickness ;
02992       ASSERT (col_thickness > 0) ;
02993       Col [col].shared1.thickness = col_thickness ;
02994
02995       /* === Remove column from degree list =========================== */
02996
02997       /* only columns in current_set will be in degree list */
02998       if (CMEMBER (col) == current_set)
02999       {
03000 #ifndef NDEBUG
03001     cols_thickness += col_thickness ;
03002 #endif
03003     cur_score = Col [col].shared2.score ;
03004     prev_col = Col [col].shared3.prev ;
03005     next_col = Col [col].shared4.degree_next ;
03006     DEBUG3 (("        cur_score "ID" prev_col "ID" next_col "ID"\n",
03007       cur_score, prev_col, next_col)) ;
03008     ASSERT (cur_score >= 0) ;
03009     ASSERT (cur_score <= n_col) ;
03010     ASSERT (cur_score >= EMPTY) ;
03011     if (prev_col == EMPTY)
03012     {
03013         head [cur_score] = next_col ;
03014     }
03015     else
03016     {
03017         Col [prev_col].shared4.degree_next = next_col ;
03018     }
03019     if (next_col != EMPTY)
03020     {
03021         Col [next_col].shared3.prev = prev_col ;
03022     }
03023       }
03024
03025       /* === Scan the column ========================================== */
03026
03027       cp = &A [Col [col].start] ;
03028       cp_end = cp + Col [col].length ;
03029       while (cp < cp_end)
03030       {
03031     /* get a row */
03032     row = *cp++ ;
03033     row_mark = Row [row].shared2.mark ;
03034     /* skip if dead */
03035     if (ROW_IS_MARKED_DEAD (row_mark))
03036     {
03037         continue ;
03038     }
03039     ASSERT (row != pivot_row) ;
03040     set_difference = row_mark - tag_mark ;
03041     /* check if the row has been seen yet */
03042     if (set_difference < 0)
03043     {
03044         ASSERT (Row [row].shared1.degree <= max_deg) ;
03045         set_difference = Row [row].shared1.degree ;
03046     }
03047     /* subtract column thickness from this row's set difference */
03048     set_difference -= col_thickness ;
03049     ASSERT (set_difference >= 0) ;
03050     /* absorb this row if the set difference becomes zero */
03051     if (set_difference == 0 && aggressive)
03052     {
03053         DEBUG3 (("aggressive absorption. Row: "ID"\n", row)) ;
03054
03055                     if (Row [row].front != EMPTY)
03056                     {
03057                         /* Row [row].front is a child of current front. */
03058                         child = Row [row].front ;
03059                         Front_parent [child] = nfr ;
03060                         DEBUG1 (("Front "ID" => front "ID", aggressive\n",
03061                                     child, nfr)) ;
03062                     }
03063                     else
03064                     {
03065                         /* this is an original row.  Keep track of which front
03066                          * assembles it, for the row-merge tree */
03067                         InFront [row] = nfr ;
03068                         DEBUG1 (("Row "ID" => front "ID", aggressive\n",
03069                                     row, nfr)) ;
03070                     }
03071
03072                     KILL_ROW (row) ;
03073
03074                     /* sum the thicknesses of all the rows */
03075                     pivot_row_thickness += Row [row].thickness ;
03076                     Row [row].thickness = 0 ;
03077     }
03078     else
03079     {
03080         /* save the new mark */
03081         Row [row].shared2.mark = set_difference + tag_mark ;
03082     }
03083       }
03084   }
03085
03086 #ifndef NDEBUG
03087   debug_deg_lists (n_row, n_col, Row, Col, head, min_score,
03088   cset_start [current_set+1]-(k+deadcol)-(cols_thickness),
03089     max_deg) ;
03090   cols_thickness = 0 ;
03091 #endif
03092
03093   /* === Add up set differences for each column ======================= */
03094
03095   DEBUG3 (("** Adding set differences phase. **\n")) ;
03096
03097   /* for each column in pivot row */
03098   rp = &A [pivot_row_start] ;
03099   rp_end = rp + pivot_row_length ;
03100   while (rp < rp_end)
03101   {
03102       /* get a column */
03103       col = *rp++ ;
03104       ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ;
03105       hash = 0 ;
03106       cur_score = 0 ;
03107       cp = &A [Col [col].start] ;
03108       /* compact the column */
03109       new_cp = cp ;
03110       cp_end = cp + Col [col].length ;
03111
03112       DEBUG4 (("Adding set diffs for Col: "ID".\n", col)) ;
03113
03114       while (cp < cp_end)
03115       {
03116     /* get a row */
03117     row = *cp++ ;
03118     ASSERT (row >= 0 && row < n_row) ;
03119     row_mark = Row [row].shared2.mark ;
03120     /* skip if dead */
03121     if (ROW_IS_MARKED_DEAD (row_mark))
03122     {
03123         DEBUG4 ((" Row "ID", dead\n", row)) ;
03124         continue ;
03125     }
03126     DEBUG4 ((" Row "ID", set diff "ID"\n", row, row_mark-tag_mark));
03127                 ASSERT (row_mark >= tag_mark) ;
03128     /* compact the column */
03129     *new_cp++ = row ;
03130     /* compute hash function */
03131     hash += row ;
03132     /* add set difference */
03133     cur_score += row_mark - tag_mark ;
03134     /* integer overflow... */
03135     cur_score = MIN (cur_score, n_col) ;
03136       }
03137
03138       /* recompute the column's length */
03139       Col [col].length = (Int) (new_cp - &A [Col [col].start]) ;
03140
03141       /* === Further mass elimination ================================= */
03142
03143       if (Col [col].length == 0 && CMEMBER (col) == current_set)
03144       {
03145     DEBUG4 (("further mass elimination. Col: "ID"\n", col)) ;
03146     /* nothing left but the pivot row in this column */
03147     KILL_PRINCIPAL_COL (col) ;
03148     pivot_row_degree -= Col [col].shared1.thickness ;
03149     ASSERT (pivot_row_degree >= 0) ;
03150     /* order it */
03151     Col [col].shared2.order = k ;
03152     /* increment order count by column thickness */
03153     k += Col [col].shared1.thickness ;
03154                 pivot_col_thickness += Col [col].shared1.thickness ;
03155                 /* add to column list of front */
03156 #ifndef NDEBUG
03157                 DEBUG1 (("Mass")) ;
03158                 dump_super (col, Col, n_col) ;
03159 #endif
03160                 Col [Col [col].lastcol].nextcol = Front_cols [nfr] ;
03161                 Front_cols [nfr] = col ;
03162       }
03163       else
03164       {
03165     /* === Prepare for supercolumn detection ==================== */
03166
03167     DEBUG4 (("Preparing supercol detection for Col: "ID".\n", col));
03168
03169     /* save score so far */
03170     Col [col].shared2.score = cur_score ;
03171
03172     /* add column to hash table, for supercolumn detection */
03173     hash %= n_col + 1 ;
03174
03175     DEBUG4 ((" Hash = "ID", n_col = "ID".\n", hash, n_col)) ;
03176     ASSERT (((Int) hash) <= n_col) ;
03177
03178     head_column = head [hash] ;
03179     if (head_column > EMPTY)
03180     {
03181         /* degree list "hash" is non-empty, use prev (shared3) of */
03182         /* first column in degree list as head of hash bucket */
03183         first_col = Col [head_column].shared3.headhash ;
03184         Col [head_column].shared3.headhash = col ;
03185     }
03186     else
03187     {
03188         /* degree list "hash" is empty, use head as hash bucket */
03189         first_col = - (head_column + 2) ;
03190         head [hash] = - (col + 2) ;
03191     }
03192     Col [col].shared4.hash_next = first_col ;
03193
03194     /* save hash function in Col [col].shared3.hash */
03195     Col [col].shared3.hash = (Int) hash ;
03196     ASSERT (COL_IS_ALIVE (col)) ;
03197       }
03198   }
03199
03200   /* The approximate external column degree is now computed.  */
03201
03202   /* === Supercolumn detection ======================================== */
03203
03204   DEBUG3 (("** Supercolumn detection phase. **\n")) ;
03205
03206   detect_super_cols (
03207 #ifndef NDEBUG
03208     n_col, Row,
03209 #endif
03210     Col, A, head, pivot_row_start, pivot_row_length, cmember) ;
03211
03212   /* === Kill the pivotal column ====================================== */
03213
03214   DEBUG1 ((" KILLING column detect supercols "ID" \n", pivot_col)) ;
03215   KILL_PRINCIPAL_COL (pivot_col) ;
03216
03217   /* add columns to column list of front */
03218 #ifndef NDEBUG
03219   DEBUG1 (("Pivot")) ;
03220   dump_super (pivot_col, Col, n_col) ;
03221 #endif
03222   Col [Col [pivot_col].lastcol].nextcol = Front_cols [nfr] ;
03223   Front_cols [nfr] = pivot_col ;
03224
03225   /* === Clear mark =================================================== */
03226
03227   tag_mark = clear_mark (tag_mark+max_deg+1, max_mark, n_row, Row) ;
03228
03229 #ifndef NDEBUG
03230   DEBUG3 (("check3\n")) ;
03231   debug_mark (n_row, Row, tag_mark, max_mark) ;
03232 #endif
03233
03234   /* === Finalize the new pivot row, and column scores ================ */
03235
03236   DEBUG3 (("** Finalize scores phase. **\n")) ;
03237
03238   /* for each column in pivot row */
03239   rp = &A [pivot_row_start] ;
03240   /* compact the pivot row */
03241   new_rp = rp ;
03242   rp_end = rp + pivot_row_length ;
03243   while (rp < rp_end)
03244   {
03245       col = *rp++ ;
03246       /* skip dead columns */
03247       if (COL_IS_DEAD (col))
03248       {
03249     continue ;
03250       }
03251       *new_rp++ = col ;
03252       /* add new pivot row to column */
03253       A [Col [col].start + (Col [col].length++)] = pivot_row ;
03254
03255       /* retrieve score so far and add on pivot row's degree. */
03256       /* (we wait until here for this in case the pivot */
03257       /* row's degree was reduced due to mass elimination). */
03258       cur_score = Col [col].shared2.score + pivot_row_degree ;
03259
03260       /* calculate the max possible score as the number of */
03261       /* external columns minus the 'k' value minus the */
03262       /* columns thickness */
03263       max_score = n_col - k - Col [col].shared1.thickness ;
03264
03265       /* make the score the external degree of the union-of-rows */
03266       cur_score -= Col [col].shared1.thickness ;
03267
03268       /* make sure score is less or equal than the max score */
03269       cur_score = MIN (cur_score, max_score) ;
03270       ASSERT (cur_score >= 0) ;
03271
03272       /* store updated score */
03273       Col [col].shared2.score = cur_score ;
03274
03275       /* === Place column back in degree list ========================= */
03276
03277       if (CMEMBER (col) == current_set)
03278       {
03279     ASSERT (min_score >= 0) ;
03280     ASSERT (min_score <= n_col) ;
03281     ASSERT (cur_score >= 0) ;
03282     ASSERT (cur_score <= n_col) ;
03283     ASSERT (head [cur_score] >= EMPTY) ;
03284     next_col = head [cur_score] ;
03285     Col [col].shared4.degree_next = next_col ;
03286     Col [col].shared3.prev = EMPTY ;
03287     if (next_col != EMPTY)
03288     {
03289         Col [next_col].shared3.prev = col ;
03290     }
03291     head [cur_score] = col ;
03292     /* see if this score is less than current min */
03293     min_score = MIN (min_score, cur_score) ;
03294       }
03295       else
03296       {
03297     Col [col].shared4.degree_next = EMPTY ;
03298     Col [col].shared3.prev = EMPTY ;
03299       }
03300   }
03301
03302 #ifndef NDEBUG
03303   debug_deg_lists (n_row, n_col, Row, Col, head,
03304     min_score, cset_start [current_set+1]-(k+deadcol), max_deg) ;
03305 #endif
03306
03307   /* frontal matrix can have more pivot cols than pivot rows for */
03308   /* singular matrices. */
03309
03310   /* number of candidate pivot columns */
03311   Front_npivcol [nfr] = pivot_col_thickness ;
03312
03313   /* all rows (not just size of contrib. block) */
03314   Front_nrows [nfr] = pivot_row_thickness ;
03315
03316   /* all cols */
03317   Front_ncols [nfr] = pivot_col_thickness + pivot_row_degree ;
03318
03319   Front_parent [nfr] = EMPTY ;
03320
03321   pivot_row_thickness -= pivot_col_thickness ;
03322   DEBUG1 (("Front "ID" Pivot_row_thickness after pivot cols elim: "ID"\n",
03323        nfr, pivot_row_thickness)) ;
03324   pivot_row_thickness = MAX (0, pivot_row_thickness) ;
03325
03326   /* === Resurrect the new pivot row ================================== */
03327
03328   if ((pivot_row_degree > 0 && pivot_row_thickness > 0 && (order_for_lu))
03329      || (pivot_row_degree > 0 && (!order_for_lu)))
03330   {
03331       /* update pivot row length to reflect any cols that were killed */
03332       /* during super-col detection and mass elimination */
03333       Row [pivot_row].start  = pivot_row_start ;
03334       Row [pivot_row].length = (Int) (new_rp - &A[pivot_row_start]) ;
03335       Row [pivot_row].shared1.degree = pivot_row_degree ;
03336       Row [pivot_row].shared2.mark = 0 ;
03337       Row [pivot_row].thickness = pivot_row_thickness ;
03338       Row [pivot_row].front = nfr ;
03339       /* pivot row is no longer dead */
03340       DEBUG1 (("Resurrect Pivot_row "ID" deg: "ID"\n",
03341       pivot_row, pivot_row_degree)) ;
03342   }
03343
03344 #ifndef NDEBUG
03345   DEBUG1 (("Front "ID" : "ID" "ID" "ID" ", nfr,
03346      Front_npivcol [nfr], Front_nrows [nfr], Front_ncols [nfr])) ;
03347   DEBUG1 ((" cols:[ ")) ;
03348   debug_d = 0 ;
03349   for (col = Front_cols [nfr] ; col != EMPTY ; col = Col [col].nextcol)
03350   {
03351     DEBUG1 ((" "ID, col)) ;
03352     ASSERT (col >= 0 && col < n_col) ;
03353     ASSERT (COL_IS_DEAD (col)) ;
03354     debug_d++ ;
03355     ASSERT (debug_d <= pivot_col_thickness) ;
03356   }
03357   ASSERT (debug_d == pivot_col_thickness) ;
03358   DEBUG1 ((" ]\n ")) ;
03359 #endif
03360    nfr++ ; /* one more front */
03361     }
03362
03363     /* === All principal columns have now been ordered ====================== */
03364
03365     *p_nfr = nfr ;
03366     return (ngarbage) ;
03367 }
03368
03369
03370 /* ========================================================================== */
03371 /* === detect_super_cols ==================================================== */
03372 /* ========================================================================== */
03373
03374 /*
03375  *  Detects supercolumns by finding matches between columns in the hash buckets.
03376  *  Check amongst columns in the set A [row_start ... row_start + row_length-1].
03377  *  The columns under consideration are currently *not* in the degree lists,
03378  *  and have already been placed in the hash buckets.
03379  *
03380  *  The hash bucket for columns whose hash function is equal to h is stored
03381  *  as follows:
03382  *
03383  *  if head [h] is >= 0, then head [h] contains a degree list, so:
03384  *
03385  *    head [h] is the first column in degree bucket h.
03386  *    Col [head [h]].headhash gives the first column in hash bucket h.
03387  *
03388  *  otherwise, the degree list is empty, and:
03389  *
03390  *    -(head [h] + 2) is the first column in hash bucket h.
03391  *
03392  *  For a column c in a hash bucket, Col [c].shared3.prev is NOT a "previous
03393  *  column" pointer.  Col [c].shared3.hash is used instead as the hash number
03394  *  for that column.  The value of Col [c].shared4.hash_next is the next column
03395  *  in the same hash bucket.
03396  *
03397  *  Assuming no, or "few" hash collisions, the time taken by this routine is
03398  *  linear in the sum of the sizes (lengths) of each column whose score has
03399  *  just been computed in the approximate degree computation.
03400  *  Not user-callable.
03401  */
03402
03403 PRIVATE void detect_super_cols
03404 (
03405     /* === Parameters ======================================================= */
03406
03407 #ifndef NDEBUG
03408     /* these two parameters are only needed when debugging is enabled: */
03409     Int n_col,      /* number of columns of A */
03410     CColamd_Row Row [ ],  /* of size n_row+1 */
03411 #endif
03412
03413     CColamd_Col Col [ ],  /* of size n_col+1 */
03414     Int A [ ],      /* row indices of A */
03415     Int head [ ],   /* head of degree lists and hash buckets */
03416     Int row_start,    /* pointer to set of columns to check */
03417     Int row_length,   /* number of columns to check */
03418     Int cmember [ ]   /* col -> cset mapping */
03419 )
03420 {
03421     /* === Local variables ================================================== */
03422
03423     Int hash ;      /* hash value for a column */
03424     Int *rp ;     /* pointer to a row */
03425     Int c ;     /* a column index */
03426     Int super_c ;   /* column index of the column to absorb into */
03427     Int *cp1 ;      /* column pointer for column super_c */
03428     Int *cp2 ;      /* column pointer for column c */
03429     Int length ;    /* length of column super_c */
03430     Int prev_c ;    /* column preceding c in hash bucket */
03431     Int i ;     /* loop counter */
03432     Int *rp_end ;   /* pointer to the end of the row */
03433     Int col ;     /* a column index in the row to check */
03434     Int head_column ;   /* first column in hash bucket or degree list */
03435     Int first_col ;   /* first column in hash bucket */
03436
03437     /* === Consider each column in the row ================================== */
03438
03439     rp = &A [row_start] ;
03440     rp_end = rp + row_length ;
03441     while (rp < rp_end)
03442     {
03443   col = *rp++ ;
03444   if (COL_IS_DEAD (col))
03445   {
03446       continue ;
03447   }
03448
03449   /* get hash number for this column */
03450   hash = Col [col].shared3.hash ;
03451   ASSERT (hash <= n_col) ;
03452
03453   /* === Get the first column in this hash bucket ===================== */
03454
03455   head_column = head [hash] ;
03456   if (head_column > EMPTY)
03457   {
03458       first_col = Col [head_column].shared3.headhash ;
03459   }
03460   else
03461   {
03462       first_col = - (head_column + 2) ;
03463   }
03464
03465   /* === Consider each column in the hash bucket ====================== */
03466
03467   for (super_c = first_col ; super_c != EMPTY ;
03468       super_c = Col [super_c].shared4.hash_next)
03469   {
03470       ASSERT (COL_IS_ALIVE (super_c)) ;
03471       ASSERT (Col [super_c].shared3.hash == hash) ;
03472       length = Col [super_c].length ;
03473
03474       /* prev_c is the column preceding column c in the hash bucket */
03475       prev_c = super_c ;
03476
03477       /* === Compare super_c with all columns after it ================ */
03478
03479       for (c = Col [super_c].shared4.hash_next ;
03480      c != EMPTY ; c = Col [c].shared4.hash_next)
03481       {
03482     ASSERT (c != super_c) ;
03483     ASSERT (COL_IS_ALIVE (c)) ;
03484     ASSERT (Col [c].shared3.hash == hash) ;
03485
03486     /* not identical if lengths or scores are different, */
03487     /* or if in different constraint sets */
03488     if (Col [c].length != length ||
03489         Col [c].shared2.score != Col [super_c].shared2.score
03490         || CMEMBER (c) != CMEMBER (super_c))
03491     {
03492         prev_c = c ;
03493         continue ;
03494     }
03495
03496     /* compare the two columns */
03497     cp1 = &A [Col [super_c].start] ;
03498     cp2 = &A [Col [c].start] ;
03499
03500     for (i = 0 ; i < length ; i++)
03501     {
03502         /* the columns are "clean" (no dead rows) */
03503         ASSERT (ROW_IS_ALIVE (*cp1)) ;
03504         ASSERT (ROW_IS_ALIVE (*cp2)) ;
03505         /* row indices will same order for both supercols, */
03506         /* no gather scatter nessasary */
03507         if (*cp1++ != *cp2++)
03508         {
03509       break ;
03510         }
03511     }
03512
03513     /* the two columns are different if the for-loop "broke" */
03514           /* super columns should belong to the same constraint set */
03515     if (i != length)
03516     {
03517         prev_c = c ;
03518         continue ;
03519     }
03520
03521     /* === Got it!  two columns are identical =================== */
03522
03523     ASSERT (Col [c].shared2.score == Col [super_c].shared2.score) ;
03524
03525     Col [super_c].shared1.thickness += Col [c].shared1.thickness ;
03526     Col [c].shared1.parent = super_c ;
03527     KILL_NON_PRINCIPAL_COL (c) ;
03528     /* order c later, in order_children() */
03529     Col [c].shared2.order = EMPTY ;
03530     /* remove c from hash bucket */
03531     Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ;
03532
03533     /* add c to end of list of super_c */
03534     ASSERT (Col [super_c].lastcol >= 0) ;
03535     ASSERT (Col [super_c].lastcol < n_col) ;
03536     Col [Col [super_c].lastcol].nextcol = c ;
03537     Col [super_c].lastcol = Col [c].lastcol ;
03538 #ifndef NDEBUG
03539     /* dump the supercolumn */
03540     DEBUG1 (("Super")) ;
03541     dump_super (super_c, Col, n_col) ;
03542 #endif
03543       }
03544   }
03545
03546   /* === Empty this hash bucket ======================================= */
03547
03548   if (head_column > EMPTY)
03549   {
03550       /* corresponding degree list "hash" is not empty */
03551       Col [head_column].shared3.headhash = EMPTY ;
03552   }
03553   else
03554   {
03555       /* corresponding degree list "hash" is empty */
03556       head [hash] = EMPTY ;
03557   }
03558     }
03559 }
03560
03561
03562 /* ========================================================================== */
03563 /* === garbage_collection =================================================== */
03564 /* ========================================================================== */
03565
03566 /*
03567  *  Defragments and compacts columns and rows in the workspace A.  Used when
03568  *  all avaliable memory has been used while performing row merging.  Returns
03569  *  the index of the first free position in A, after garbage collection.  The
03570  *  time taken by this routine is linear is the size of the array A, which is
03571  *  itself linear in the number of nonzeros in the input matrix.
03572  *  Not user-callable.
03573  */
03574
03575 PRIVATE Int garbage_collection  /* returns the new value of pfree */
03576 (
03577     /* === Parameters ======================================================= */
03578
03579     Int n_row,      /* number of rows */
03580     Int n_col,      /* number of columns */
03581     CColamd_Row Row [ ],  /* row info */
03582     CColamd_Col Col [ ],  /* column info */
03583     Int A [ ],      /* A [0 ... Alen-1] holds the matrix */
03584     Int *pfree      /* &A [0] ... pfree is in use */
03585 )
03586 {
03587     /* === Local variables ================================================== */
03588
03589     Int *psrc ;     /* source pointer */
03590     Int *pdest ;    /* destination pointer */
03591     Int j ;     /* counter */
03592     Int r ;     /* a row index */
03593     Int c ;     /* a column index */
03594     Int length ;    /* length of a row or column */
03595
03596 #ifndef NDEBUG
03597     Int debug_rows ;
03598     DEBUG2 (("Defrag..\n")) ;
03599     for (psrc = &A[0] ; psrc < pfree ; psrc++) ASSERT (*psrc >= 0) ;
03600     debug_rows = 0 ;
03601 #endif
03602
03603     /* === Defragment the columns =========================================== */
03604
03605     pdest = &A[0] ;
03606     for (c = 0 ; c < n_col ; c++)
03607     {
03608   if (COL_IS_ALIVE (c))
03609   {
03610       psrc = &A [Col [c].start] ;
03611
03612       /* move and compact the column */
03613       ASSERT (pdest <= psrc) ;
03614       Col [c].start = (Int) (pdest - &A [0]) ;
03615       length = Col [c].length ;
03616       for (j = 0 ; j < length ; j++)
03617       {
03618     r = *psrc++ ;
03619     if (ROW_IS_ALIVE (r))
03620     {
03621         *pdest++ = r ;
03622     }
03623       }
03624       Col [c].length = (Int) (pdest - &A [Col [c].start]) ;
03625   }
03626     }
03627
03628     /* === Prepare to defragment the rows =================================== */
03629
03630     for (r = 0 ; r < n_row ; r++)
03631     {
03632   if (ROW_IS_DEAD (r) || (Row [r].length == 0))
03633   {
03634       /* This row is already dead, or is of zero length.  Cannot compact
03635        * a row of zero length, so kill it.  NOTE: in the current version,
03636        * there are no zero-length live rows.  Kill the row (for the first
03637        * time, or again) just to be safe. */
03638       KILL_ROW (r) ;
03639   }
03640   else
03641   {
03642       /* save first column index in Row [r].shared2.first_column */
03643       psrc = &A [Row [r].start] ;
03644       Row [r].shared2.first_column = *psrc ;
03645       ASSERT (ROW_IS_ALIVE (r)) ;
03646       /* flag the start of the row with the one's complement of row */
03647       *psrc = ONES_COMPLEMENT (r) ;
03648 #ifndef NDEBUG
03649       debug_rows++ ;
03650 #endif
03651   }
03652     }
03653
03654     /* === Defragment the rows ============================================== */
03655
03656     psrc = pdest ;
03657     while (psrc < pfree)
03658     {
03659   /* find a negative number ... the start of a row */
03660   if (*psrc++ < 0)
03661   {
03662       psrc-- ;
03663       /* get the row index */
03664       r = ONES_COMPLEMENT (*psrc) ;
03665       ASSERT (r >= 0 && r < n_row) ;
03666       /* restore first column index */
03667       *psrc = Row [r].shared2.first_column ;
03668       ASSERT (ROW_IS_ALIVE (r)) ;
03669
03670       /* move and compact the row */
03671       ASSERT (pdest <= psrc) ;
03672       Row [r].start = (Int) (pdest - &A [0]) ;
03673       length = Row [r].length ;
03674       for (j = 0 ; j < length ; j++)
03675       {
03676     c = *psrc++ ;
03677     if (COL_IS_ALIVE (c))
03678     {
03679         *pdest++ = c ;
03680     }
03681       }
03682       Row [r].length = (Int) (pdest - &A [Row [r].start]) ;
03683 #ifndef NDEBUG
03684       debug_rows-- ;
03685 #endif
03686   }
03687     }
03688
03689     /* ensure we found all the rows */
03690     ASSERT (debug_rows == 0) ;
03691
03692     /* === Return the new value of pfree ==================================== */
03693
03694     return ((Int) (pdest - &A [0])) ;
03695 }
03696
03697
03698 /* ========================================================================== */
03699 /* === clear_mark =========================================================== */
03700 /* ========================================================================== */
03701
03702 /*
03703  *  Clears the Row [ ].shared2.mark array, and returns the new tag_mark.
03704  *  Return value is the new tag_mark.  Not user-callable.
03705  */
03706
03707 PRIVATE Int clear_mark  /* return the new value for tag_mark */
03708 (
03709     /* === Parameters ======================================================= */
03710
03711     Int tag_mark, /* new value of tag_mark */
03712     Int max_mark, /* max allowed value of tag_mark */
03713
03714     Int n_row,    /* number of rows in A */
03715     CColamd_Row Row [ ] /* Row [0 ... n_row-1].shared2.mark is set to zero */
03716 )
03717 {
03718     /* === Local variables ================================================== */
03719
03720     Int r ;
03721
03722     if (tag_mark <= 0 || tag_mark >= max_mark)
03723     {
03724   for (r = 0 ; r < n_row ; r++)
03725   {
03726       if (ROW_IS_ALIVE (r))
03727       {
03728     Row [r].shared2.mark = 0 ;
03729       }
03730   }
03731   tag_mark = 1 ;
03732     }
03733
03734     return (tag_mark) ;
03735 }
03736
03737
03738 /* ========================================================================== */
03739 /* === print_report ========================================================= */
03740 /* ========================================================================== */
03741
03742 /* No printing occurs if NPRINT is defined at compile time. */
03743
03744 PRIVATE void print_report
03745 (
03746     char *method,
03747     Int stats [CCOLAMD_STATS]
03748 )
03749 {
03750
03751     Int i1, i2, i3 ;
03752
03753     PRINTF (("\n%s version %d.%d, %s: ", method,
03754       CCOLAMD_MAIN_VERSION, CCOLAMD_SUB_VERSION, CCOLAMD_DATE)) ;
03755
03756     if (!stats)
03757     {
03758       PRINTF (("No statistics available.\n")) ;
03759   return ;
03760     }
03761
03762     i1 = stats [CCOLAMD_INFO1] ;
03763     i2 = stats [CCOLAMD_INFO2] ;
03764     i3 = stats [CCOLAMD_INFO3] ;
03765
03766     if (stats [CCOLAMD_STATUS] >= 0)
03767     {
03768       PRINTF(("OK.  ")) ;
03769     }
03770     else
03771     {
03772       PRINTF(("ERROR.  ")) ;
03773     }
03774
03775     switch (stats [CCOLAMD_STATUS])
03776     {
03777
03778   case CCOLAMD_OK_BUT_JUMBLED:
03779
03780       PRINTF(("Matrix has unsorted or duplicate row indices.\n")) ;
03781
03782       PRINTF(("%s: duplicate or out-of-order row indices:    "ID"\n",
03783         method, i3)) ;
03784
03785       PRINTF(("%s: last seen duplicate or out-of-order row:  "ID"\n",
03786         method, INDEX (i2))) ;
03787
03788       PRINTF(("%s: last seen in column:                      "ID"",
03789         method, INDEX (i1))) ;
03790
03791       /* no break - fall through to next case instead */
03792
03793   case CCOLAMD_OK:
03794
03795       PRINTF(("\n")) ;
03796
03797       PRINTF(("%s: number of dense or empty rows ignored:    "ID"\n",
03798         method, stats [CCOLAMD_DENSE_ROW])) ;
03799
03800       PRINTF(("%s: number of dense or empty columns ignored: "ID"\n",
03801         method, stats [CCOLAMD_DENSE_COL])) ;
03802
03803       PRINTF(("%s: number of garbage collections performed:  "ID"\n",
03804         method, stats [CCOLAMD_DEFRAG_COUNT])) ;
03805       break ;
03806
03807   case CCOLAMD_ERROR_A_not_present:
03808
03809       PRINTF(("Array A (row indices of matrix) not present.\n")) ;
03810       break ;
03811
03812   case CCOLAMD_ERROR_p_not_present:
03813
03814       PRINTF(("Array p (column pointers for matrix) not present.\n")) ;
03815       break ;
03816
03817   case CCOLAMD_ERROR_nrow_negative:
03818
03819       PRINTF(("Invalid number of rows ("ID").\n", i1)) ;
03820       break ;
03821
03822   case CCOLAMD_ERROR_ncol_negative:
03823
03824       PRINTF(("Invalid number of columns ("ID").\n", i1)) ;
03825       break ;
03826
03827   case CCOLAMD_ERROR_nnz_negative:
03828
03829       PRINTF(("Invalid number of nonzero entries ("ID").\n", i1)) ;
03830       break ;
03831
03832   case CCOLAMD_ERROR_p0_nonzero:
03833
03834       PRINTF(("Invalid column pointer, p [0] = "ID", must be 0.\n", i1)) ;
03835       break ;
03836
03837   case CCOLAMD_ERROR_A_too_small:
03838
03839       PRINTF(("Array A too small.\n")) ;
03840       PRINTF(("        Need Alen >= "ID", but given only Alen = "ID".\n",
03841         i1, i2)) ;
03842       break ;
03843
03844   case CCOLAMD_ERROR_col_length_negative:
03845
03846       PRINTF(("Column "ID" has a negative number of entries ("ID").\n",
03847         INDEX (i1), i2)) ;
03848       break ;
03849
03850   case CCOLAMD_ERROR_row_index_out_of_bounds:
03851
03852       PRINTF(("Row index (row "ID") out of bounds ("ID" to "ID") in"
03853         "column "ID".\n", INDEX (i2), INDEX (0), INDEX (i3-1),
03854         INDEX (i1))) ;
03855       break ;
03856
03857   case CCOLAMD_ERROR_out_of_memory:
03858
03859       PRINTF(("Out of memory.\n")) ;
03860       break ;
03861
03862   case CCOLAMD_ERROR_invalid_cmember:
03863
03864       PRINTF(("cmember invalid\n")) ;
03865       break ;
03866     }
03867 }
03868
03869
03870 /* ========================================================================= */
03871 /* === "Expert" routines =================================================== */
03872 /* ========================================================================= */
03873
03874 /* The following routines are visible outside this routine, but are not meant
03875  * to be called by the user.  They are meant for a future version of UMFPACK,
03876  * to replace UMFPACK internal routines with a similar name.
03877  */
03878
03879
03880 /* ========================================================================== */
03881 /* === CCOLAMD_apply_order ================================================== */
03882 /* ========================================================================== */
03883
03884 /*
03885  * Apply post-ordering of supernodal elimination tree.
03886  */
03887
03888 GLOBAL void CCOLAMD_apply_order
03889 (
03890     Int Front [ ],      /* of size nn on input, size nfr on output */
03891     const Int Order [ ],    /* Order [i] = k, i in the range 0..nn-1,
03892            * and k in the range 0..nfr-1, means that node
03893            * i is the kth node in the postordered tree. */
03894     Int Temp [ ],     /* workspace of size nfr */
03895     Int nn,       /* nodes are numbered in the range 0..nn-1 */
03896     Int nfr       /* the number of nodes actually in use */
03897 )
03898 {
03899     Int i, k ;
03900     for (i = 0 ; i < nn ; i++)
03901     {
03902   k = Order [i] ;
03903   ASSERT (k >= EMPTY && k < nfr) ;
03904   if (k != EMPTY)
03905   {
03906       Temp [k] = Front [i] ;
03907   }
03908     }
03909
03910     for (k = 0 ; k < nfr ; k++)
03911     {
03912   Front [k] = Temp [k] ;
03913     }
03914 }
03915
03916
03917 /* ========================================================================== */
03918 /* === CCOLAMD_fsize ======================================================== */
03919 /* ========================================================================== */
03920
03921 /* Determine the largest frontal matrix size for each subtree.
03922  * Only required to sort the children of each
03923  * node prior to postordering the column elimination tree. */
03924
03925 GLOBAL void CCOLAMD_fsize
03926 (
03927     Int nn,
03928     Int Fsize [ ],
03929     Int Fnrows [ ],
03930     Int Fncols [ ],
03931     Int Parent [ ],
03932     Int Npiv [ ]
03933 )
03934 {
03935     double dr, dc ;
03936     Int j, parent, frsize, r, c ;
03937
03938     for (j = 0 ; j < nn ; j++)
03939     {
03940   Fsize [j] = EMPTY ;
03941     }
03942
03943     /* ---------------------------------------------------------------------- */
03944     /* find max front size for tree rooted at node j, for each front j */
03945     /* ---------------------------------------------------------------------- */
03946
03947     DEBUG1 (("\n\n========================================FRONTS:\n")) ;
03948     for (j = 0 ; j < nn ; j++)
03949     {
03950   if (Npiv [j] > 0)
03951   {
03952       /* this is a frontal matrix */
03953       parent = Parent [j] ;
03954       r = Fnrows [j] ;
03955       c = Fncols [j] ;
03956       /* avoid integer overflow */
03957       dr = (double) r ;
03958       dc = (double) c ;
03959       frsize = (INT_OVERFLOW (dr * dc)) ?  Int_MAX : (r * c) ;
03960       DEBUG1 ((""ID" : npiv "ID" size "ID" parent "ID" ",
03961     j, Npiv [j], frsize, parent)) ;
03962       Fsize [j] = MAX (Fsize [j], frsize) ;
03963       DEBUG1 (("Fsize [j = "ID"] = "ID"\n", j, Fsize [j])) ;
03964       if (parent != EMPTY)
03965       {
03966     /* find the maximum frontsize of self and children */
03967     ASSERT (Npiv [parent] > 0) ;
03968     ASSERT (parent > j) ;
03969     Fsize [parent] = MAX (Fsize [parent], Fsize [j]) ;
03970     DEBUG1 (("Fsize [parent = "ID"] = "ID"\n",
03971         parent, Fsize [parent]));
03972       }
03973   }
03974     }
03975     DEBUG1 (("fsize done\n")) ;
03976 }
03977
03978
03979 /* ========================================================================= */
03980 /* === CCOLAMD_postorder =================================================== */
03981 /* ========================================================================= */
03982
03983 /* Perform a postordering (via depth-first search) of an assembly tree. */
03984
03985 GLOBAL void CCOLAMD_postorder
03986 (
03987     /* inputs, not modified on output: */
03988     Int nn,   /* nodes are in the range 0..nn-1 */
03989     Int Parent [ ], /* Parent [j] is the parent of j, or EMPTY if root */
03990     Int Nv [ ],   /* Nv [j] > 0 number of pivots represented by node j,
03991        * or zero if j is not a node. */
03992     Int Fsize [ ],  /* Fsize [j]: size of node j */
03993
03994     /* output, not defined on input: */
03995     Int Order [ ],  /* output post-order */
03996
03997     /* workspaces of size nn: */
03998     Int Child [ ],
03999     Int Sibling [ ],
04000     Int Stack [ ],
04001     Int Front_cols [ ],
04002
04003     /* input, not modified on output: */
04004     Int cmember [ ]
04005 )
04006 {
04007     Int i, j, k, parent, frsize, f, fprev, maxfrsize, bigfprev, bigf, fnext ;
04008
04009     for (j = 0 ; j < nn ; j++)
04010     {
04011   Child [j] = EMPTY ;
04012   Sibling [j] = EMPTY ;
04013     }
04014
04015     /* --------------------------------------------------------------------- */
04016     /* place the children in link lists - bigger elements tend to be last */
04017     /* --------------------------------------------------------------------- */
04018
04019     for (j = nn-1 ; j >= 0 ; j--)
04020     {
04021   if (Nv [j] > 0)
04022   {
04023       /* this is an element */
04024       parent = Parent [j] ;
04025       if (parent != EMPTY)
04026       {
04027     /* place the element in link list of the children its parent */
04028     /* bigger elements will tend to be at the end of the list */
04029     Sibling [j] = Child [parent] ;
04030     if (CMEMBER (Front_cols[parent]) == CMEMBER (Front_cols[j]))
04031     {
04032         Child [parent] = j ;
04033     }
04034       }
04035   }
04036     }
04037
04038 #ifndef NDEBUG
04039     {
04040   Int nels, ff, nchild ;
04041   DEBUG1 (("\n\n================================ ccolamd_postorder:\n"));
04042   nels = 0 ;
04043   for (j = 0 ; j < nn ; j++)
04044   {
04045       if (Nv [j] > 0)
04046       {
04047     DEBUG1 ((""ID" :  nels "ID" npiv "ID" size "ID
04048         " parent "ID" maxfr "ID"\n", j, nels,
04049         Nv [j], Fsize [j], Parent [j], Fsize [j])) ;
04050     /* this is an element */
04051     /* dump the link list of children */
04052     nchild = 0 ;
04053     DEBUG1 (("    Children: ")) ;
04054     for (ff = Child [j] ; ff != EMPTY ; ff = Sibling [ff])
04055     {
04056         DEBUG1 ((ID" ", ff)) ;
04057         nchild++ ;
04058         ASSERT (nchild < nn) ;
04059     }
04060     DEBUG1 (("\n")) ;
04061     parent = Parent [j] ;
04062     nels++ ;
04063       }
04064   }
04065     }
04066 #endif
04067
04068     /* --------------------------------------------------------------------- */
04069     /* place the largest child last in the list of children for each node */
04070     /* --------------------------------------------------------------------- */
04071
04072     for (i = 0 ; i < nn ; i++)
04073     {
04074   if (Nv [i] > 0 && Child [i] != EMPTY)
04075   {
04076
04077 #ifndef NDEBUG
04078       Int nchild ;
04079       DEBUG1 (("Before partial sort, element "ID"\n", i)) ;
04080       nchild = 0 ;
04081       for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04082       {
04083     DEBUG1 (("      f: "ID"  size: "ID"\n", f, Fsize [f])) ;
04084     nchild++ ;
04085       }
04086 #endif
04087
04088       /* find the biggest element in the child list */
04089       fprev = EMPTY ;
04090       maxfrsize = EMPTY ;
04091       bigfprev = EMPTY ;
04092       bigf = EMPTY ;
04093       for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04094       {
04095     frsize = Fsize [f] ;
04096     if (frsize >= maxfrsize)
04097     {
04098         /* this is the biggest seen so far */
04099         maxfrsize = frsize ;
04100         bigfprev = fprev ;
04101         bigf = f ;
04102     }
04103     fprev = f ;
04104       }
04105
04106       fnext = Sibling [bigf] ;
04107
04108       DEBUG1 (("bigf "ID" maxfrsize "ID" bigfprev "ID" fnext "ID
04109     " fprev " ID"\n", bigf, maxfrsize, bigfprev, fnext, fprev)) ;
04110
04111       if (fnext != EMPTY)
04112       {
04113     /* if fnext is EMPTY then bigf is already at the end of list */
04114
04115     if (bigfprev == EMPTY)
04116     {
04117         /* delete bigf from the element of the list */
04118         Child [i] = fnext ;
04119     }
04120     else
04121     {
04122         /* delete bigf from the middle of the list */
04123         Sibling [bigfprev] = fnext ;
04124     }
04125
04126     /* put bigf at the end of the list */
04127     Sibling [bigf] = EMPTY ;
04128     Sibling [fprev] = bigf ;
04129       }
04130
04131 #ifndef NDEBUG
04132       DEBUG1 (("After partial sort, element "ID"\n", i)) ;
04133       for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04134       {
04135     DEBUG1 (("        "ID"  "ID"\n", f, Fsize [f])) ;
04136     nchild-- ;
04137       }
04138 #endif
04139   }
04140     }
04141
04142     /* --------------------------------------------------------------------- */
04143     /* postorder the assembly tree */
04144     /* --------------------------------------------------------------------- */
04145
04146     for (i = 0 ; i < nn ; i++)
04147     {
04148   Order [i] = EMPTY ;
04149     }
04150
04151     k = 0 ;
04152
04153     for (i = 0 ; i < nn ; i++)
04154     {
04155   if ((Parent [i] == EMPTY
04156       || (CMEMBER (Front_cols [Parent [i]]) != CMEMBER (Front_cols [i])))
04157       && Nv [i] > 0)
04158   {
04159       DEBUG1 (("Root of assembly tree "ID"\n", i)) ;
04160       k = CCOLAMD_post_tree (i, k, Child, Sibling, Order, Stack) ;
04161   }
04162     }
04163 }
04164
04165
04166 /* ========================================================================= */
04167 /* === CCOLAMD_post_tree =================================================== */
04168 /* ========================================================================= */
04169
04170 /* Post-ordering of a supernodal column elimination tree.  */
04171
04172 GLOBAL Int CCOLAMD_post_tree
04173 (
04174     Int root,     /* root of the tree */
04175     Int k,      /* start numbering at k */
04176     Int Child [ ],    /* input argument of size nn, undefined on
04177          * output.  Child [i] is the head of a link
04178          * list of all nodes that are children of node
04179          * i in the tree. */
04180     const Int Sibling [ ],  /* input argument of size nn, not modified.
04181          * If f is a node in the link list of the
04182          * children of node i, then Sibling [f] is the
04183          * next child of node i.
04184          */
04185     Int Order [ ],    /* output order, of size nn.  Order [i] = k
04186          * if node i is the kth node of the reordered
04187          * tree. */
04188     Int Stack [ ]   /* workspace of size nn */
04189 )
04190 {
04191     Int f, head, h, i ;
04192
04193 #if 0
04194     /* --------------------------------------------------------------------- */
04195     /* recursive version (Stack [ ] is not used): */
04196     /* --------------------------------------------------------------------- */
04197
04198     /* this is simple, but can cause stack overflow if nn is large */
04199     i = root ;
04200     for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04201     {
04202   k = CCOLAMD_post_tree (f, k, Child, Sibling, Order, Stack, nn) ;
04203     }
04204     Order [i] = k++ ;
04205     return (k) ;
04206 #endif
04207
04208     /* --------------------------------------------------------------------- */
04209     /* non-recursive version, using an explicit stack */
04210     /* --------------------------------------------------------------------- */
04211
04212     /* push root on the stack */
04213     head = 0 ;
04214     Stack [0] = root ;
04215
04216     while (head >= 0)
04217     {
04218   /* get head of stack */
04219   i = Stack [head] ;
04220   DEBUG1 (("head of stack "ID" \n", i)) ;
04221
04222   if (Child [i] != EMPTY)
04223   {
04224       /* the children of i are not yet ordered */
04225       /* push each child onto the stack in reverse order */
04226       /* so that small ones at the head of the list get popped first */
04227       /* and the biggest one at the end of the list gets popped last */
04228       for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04229       {
04230     head++ ;
04231       }
04232       h = head ;
04233       for (f = Child [i] ; f != EMPTY ; f = Sibling [f])
04234       {
04235     ASSERT (h > 0) ;
04236     Stack [h--] = f ;
04237     DEBUG1 (("push "ID" on stack\n", f)) ;
04238       }
04239       ASSERT (Stack [h] == i) ;
04240
04241       /* delete child list so that i gets ordered next time we see it */
04242       Child [i] = EMPTY ;
04243   }
04244   else
04245   {
04246       /* the children of i (if there were any) are already ordered */
04247       /* remove i from the stack and order it.  Front i is kth front */
04248       head-- ;
04249       DEBUG1 (("pop "ID" order "ID"\n", i, k)) ;
04250       Order [i] = k++ ;
04251   }
04252
04253 #ifndef NDEBUG
04254   DEBUG1 (("\nStack:")) ;
04255   for (h = head ; h >= 0 ; h--)
04256   {
04257       Int j = Stack [h] ;
04258       DEBUG1 ((" "ID, j)) ;
04259   }
04260   DEBUG1 (("\n\n")) ;
04261 #endif
04262
04263     }
04264     return (k) ;
04265 }
04266
04267
04268
04269 /* ========================================================================== */
04270 /* === CCOLAMD debugging routines =========================================== */
04271 /* ========================================================================== */
04272
04273 /* When debugging is disabled, the remainder of this file is ignored. */
04274
04275 #ifndef NDEBUG
04276
04277
04278 /* ========================================================================== */
04279 /* === debug_structures ===================================================== */
04280 /* ========================================================================== */
04281
04282 /*
04283  *  At this point, all empty rows and columns are dead.  All live columns
04284  *  are "clean" (containing no dead rows) and simplicial (no supercolumns
04285  *  yet).  Rows may contain dead columns, but all live rows contain at
04286  *  least one live column.
04287  */
04288
04289 PRIVATE void debug_structures
04290 (
04291     /* === Parameters ======================================================= */
04292
04293     Int n_row,
04294     Int n_col,
04295     CColamd_Row Row [ ],
04296     CColamd_Col Col [ ],
04297     Int A [ ],
04298     Int cmember [ ],
04299     Int cset_start [ ]
04300 )
04301 {
04302     /* === Local variables ================================================== */
04303
04304     Int i ;
04305     Int c ;
04306     Int *cp ;
04307     Int *cp_end ;
04308     Int len ;
04309     Int score ;
04310     Int r ;
04311     Int *rp ;
04312     Int *rp_end ;
04313     Int deg ;
04314     Int cs ;
04315
04316     /* === Check A, Row, and Col ============================================ */
04317
04318     for (c = 0 ; c < n_col ; c++)
04319     {
04320   if (COL_IS_ALIVE (c))
04321   {
04322       len = Col [c].length ;
04323       score = Col [c].shared2.score ;
04324       DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ;
04325       ASSERT (len > 0) ;
04326       ASSERT (score >= 0) ;
04327       ASSERT (Col [c].shared1.thickness == 1) ;
04328       cp = &A [Col [c].start] ;
04329       cp_end = cp + len ;
04330       while (cp < cp_end)
04331       {
04332     r = *cp++ ;
04333     ASSERT (ROW_IS_ALIVE (r)) ;
04334       }
04335   }
04336   else
04337   {
04338       i = Col [c].shared2.order ;
04339       cs = CMEMBER (c) ;
04340       ASSERT (i >= cset_start [cs] && i < cset_start [cs+1]) ;
04341   }
04342     }
04343
04344     for (r = 0 ; r < n_row ; r++)
04345     {
04346   if (ROW_IS_ALIVE (r))
04347   {
04348       i = 0 ;
04349       len = Row [r].length ;
04350       deg = Row [r].shared1.degree ;
04351       ASSERT (len > 0) ;
04352       ASSERT (deg > 0) ;
04353       rp = &A [Row [r].start] ;
04354       rp_end = rp + len ;
04355       while (rp < rp_end)
04356       {
04357     c = *rp++ ;
04358     if (COL_IS_ALIVE (c))
04359     {
04360         i++ ;
04361     }
04362       }
04363       ASSERT (i > 0) ;
04364   }
04365     }
04366 }
04367
04368
04369 /* ========================================================================== */
04370 /* === debug_deg_lists ====================================================== */
04371 /* ========================================================================== */
04372
04373 /*
04374  *  Prints the contents of the degree lists.  Counts the number of columns
04375  *  in the degree list and compares it to the total it should have.  Also
04376  *  checks the row degrees.
04377  */
04378
04379 PRIVATE void debug_deg_lists
04380 (
04381     /* === Parameters ======================================================= */
04382
04383     Int n_row,
04384     Int n_col,
04385     CColamd_Row Row [ ],
04386     CColamd_Col Col [ ],
04387     Int head [ ],
04388     Int min_score,
04389     Int should,
04390     Int max_deg
04391 )
04392
04393 {
04394     /* === Local variables ================================================== */
04395
04396     Int deg ;
04397     Int col ;
04398     Int have ;
04399     Int row ;
04400
04401     /* === Check the degree lists =========================================== */
04402
04403     if (n_col > 10000 && ccolamd_debug <= 0)
04404     {
04405   return ;
04406     }
04407     have = 0 ;
04408     DEBUG4 (("Degree lists: "ID"\n", min_score)) ;
04409     for (deg = 0 ; deg <= n_col ; deg++)
04410     {
04411   col = head [deg] ;
04412   if (col == EMPTY)
04413   {
04414       continue ;
04415   }
04416   DEBUG4 (("%d:", deg)) ;
04417   ASSERT (Col [col].shared3.prev == EMPTY) ;
04418   while (col != EMPTY)
04419   {
04420       DEBUG4 ((" "ID"", col)) ;
04421       have += Col [col].shared1.thickness ;
04422       ASSERT (COL_IS_ALIVE (col)) ;
04423       col = Col [col].shared4.degree_next ;
04424   }
04425   DEBUG4 (("\n")) ;
04426     }
04427     DEBUG4 (("should "ID" have "ID"\n", should, have)) ;
04428     ASSERT (should == have) ;
04429
04430     /* === Check the row degrees ============================================ */
04431
04432     if (n_row > 10000 && ccolamd_debug <= 0)
04433     {
04434   return ;
04435     }
04436     for (row = 0 ; row < n_row ; row++)
04437     {
04438   if (ROW_IS_ALIVE (row))
04439   {
04440       ASSERT (Row [row].shared1.degree <= max_deg) ;
04441   }
04442     }
04443 }
04444
04445
04446 /* ========================================================================== */
04447 /* === debug_mark =========================================================== */
04448 /* ========================================================================== */
04449
04450 /*
04451  *  Ensures that the tag_mark is less that the maximum and also ensures that
04452  *  each entry in the mark array is less than the tag mark.
04453  */
04454
04455 PRIVATE void debug_mark
04456 (
04457     /* === Parameters ======================================================= */
04458
04459     Int n_row,
04460     CColamd_Row Row [ ],
04461     Int tag_mark,
04462     Int max_mark
04463 )
04464 {
04465     /* === Local variables ================================================== */
04466
04467     Int r ;
04468
04469     /* === Check the Row marks ============================================== */
04470
04471     ASSERT (tag_mark > 0 && tag_mark <= max_mark) ;
04472     if (n_row > 10000 && ccolamd_debug <= 0)
04473     {
04474   return ;
04475     }
04476     for (r = 0 ; r < n_row ; r++)
04477     {
04478   ASSERT (Row [r].shared2.mark < tag_mark) ;
04479     }
04480 }
04481
04482
04483 /* ========================================================================== */
04484 /* === debug_matrix ========================================================= */
04485 /* ========================================================================== */
04486
04487 /* Prints out the contents of the columns and the rows.  */
04488
04489 PRIVATE void debug_matrix
04490 (
04491     /* === Parameters ======================================================= */
04492
04493     Int n_row,
04494     Int n_col,
04495     CColamd_Row Row [ ],
04496     CColamd_Col Col [ ],
04497     Int A [ ]
04498 )
04499 {
04500     /* === Local variables ================================================== */
04501
04502     Int r ;
04503     Int c ;
04504     Int *rp ;
04505     Int *rp_end ;
04506     Int *cp ;
04507     Int *cp_end ;
04508
04509     /* === Dump the rows and columns of the matrix ========================== */
04510
04511     if (ccolamd_debug < 3)
04512     {
04513   return ;
04514     }
04515     DEBUG3 (("DUMP MATRIX:\n")) ;
04516     for (r = 0 ; r < n_row ; r++)
04517     {
04518   DEBUG3 (("Row "ID" alive? "ID"\n", r, ROW_IS_ALIVE (r))) ;
04519   if (ROW_IS_DEAD (r))
04520   {
04521       continue ;
04522   }
04523
04524   DEBUG3 (("start "ID" length "ID" degree "ID"\nthickness "ID"\n",
04525     Row [r].start, Row [r].length, Row [r].shared1.degree,
04526     Row [r].thickness)) ;
04527
04528   rp = &A [Row [r].start] ;
04529   rp_end = rp + Row [r].length ;
04530   while (rp < rp_end)
04531   {
04532       c = *rp++ ;
04533       DEBUG4 (("  "ID" col "ID"\n", COL_IS_ALIVE (c), c)) ;
04534   }
04535     }
04536
04537     for (c = 0 ; c < n_col ; c++)
04538     {
04539   DEBUG3 (("Col "ID" alive? "ID"\n", c, COL_IS_ALIVE (c))) ;
04540   if (COL_IS_DEAD (c))
04541   {
04542       continue ;
04543   }
04544   DEBUG3 (("start "ID" length "ID" shared1 "ID" shared2 "ID"\n",
04545     Col [c].start, Col [c].length,
04546     Col [c].shared1.thickness, Col [c].shared2.score)) ;
04547   cp = &A [Col [c].start] ;
04548   cp_end = cp + Col [c].length ;
04549   while (cp < cp_end)
04550   {
04551       r = *cp++ ;
04552       DEBUG4 (("  "ID" row "ID"\n", ROW_IS_ALIVE (r), r)) ;
04553   }
04554     }
04555 }
04556
04557
04558 /* ========================================================================== */
04559 /* === dump_super =========================================================== */
04560 /* ========================================================================== */
04561
04562 PRIVATE void dump_super
04563 (
04564     Int super_c,
04565     CColamd_Col Col [ ],
04566     Int n_col
04567 )
04568 {
04569     Int col, ncols ;
04570
04571     DEBUG1 ((" =[ ")) ;
04572     ncols = 0 ;
04573     for (col = super_c ; col != EMPTY ; col = Col [col].nextcol)
04574     {
04575         DEBUG1 ((" "ID, col)) ;
04576         ASSERT (col >= 0 && col < n_col) ;
04577         if (col != super_c)
04578         {
04579             ASSERT (COL_IS_DEAD (col)) ;
04580         }
04581         if (Col [col].nextcol == EMPTY)
04582         {
04583             ASSERT (col == Col [super_c].lastcol) ;
04584         }
04585         ncols++ ;
04586         ASSERT (ncols <= Col [super_c].shared1.thickness) ;
04587     }
04588     ASSERT (ncols == Col [super_c].shared1.thickness) ;
04589     DEBUG1 (("]\n")) ;
04590 }
04591
04592
04593 /* ========================================================================== */
04594 /* === ccolamd_get_debug ==================================================== */
04595 /* ========================================================================== */
04596
04597 PRIVATE void ccolamd_get_debug
04598 (
04599     char *method
04600 )
04601 {
04602     FILE *debug_file ;
04603     ccolamd_debug = 0 ;   /* no debug printing */
04604
04605     /* Read debug info from the debug file. */
04606     debug_file = fopen ("debug", "r") ;
04607     if (debug_file)
04608     {
04609   (void) fscanf (debug_file, ""ID"", &ccolamd_debug) ;
04610   (void) fclose (debug_file) ;
04611     }
04612
04613     DEBUG0 ((":")) ;
04614     DEBUG1 (("%s: debug version, D = "ID" (THIS WILL BE SLOW!)\n",
04615       method, ccolamd_debug)) ;
04616     DEBUG1 ((" Debug printing level: "ID"\n", ccolamd_debug)) ;
04617 }
04618
04619 #endif