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